tutorials/pyMBE_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # pyMBE Tutorial: the Python-based Molecule Builder for ESPResSo.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective of pyMBE is to facilitate the set-up of complex molecules into the Molecular Dynamics software [espressomd.org](ESPResSo). pyMBE is specially well-suitedto facilitate setting up constant pH and grand-reaction ensemble simulations in ESPResSo. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents:\n",
    "* [Introduction](#introduction)\n",
    "* [How to create particles](#particles)\n",
    "* [How to create simple polymers](#simple_polymers)\n",
    "* [How to create complex polymers](#complex_polymers)\n",
    "* [How to create di-block copolymers](#diblock_copolymers)\n",
    "* [Practice by creating an alternating copolymer](#alternating_copolymers)\n",
    "* [How to create peptides](#peptides)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "\n",
    "Let us get started by importing pyMBE library and other important libraries for this tutorial, such as ESPResSo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import pyMBE and  ESPResSo\n",
    "import pyMBE\n",
    "import espressomd\n",
    "from espressomd import interactions\n",
    "\n",
    "# Only necesary to produce the pictures used in this tutorial\n",
    "from lib.handy_functions import do_snapshot_espresso_system\n",
    "from PIL import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating an instance of pyMBE we can make our code shorter and more readable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb = pyMBE.pymbe_library(seed=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When pyMBE is inicialized, a default system of reduced units is defined. \n",
    "\n",
    "* Unit_length = 0.355 nm.\n",
    "* Unit_charge = 1 elementary charge.\n",
    "* Temperature = 298.15 K.\n",
    "\n",
    "The active set of reduced units can be consulted using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_unit_set = pmb.get_reduced_units()\n",
    "print(reduced_unit_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This default definition of reduced units can be changed at the convience of the user using the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.set_reduced_units(unit_length = 0.5*pmb.units.nm,  \n",
    "                      unit_charge = 5*pmb.units.e)\n",
    "                      #, temperature=300*pmb.units.K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: All input variables will be given to ESPResSo using these reduced units, since it is a convinient choice for the simulation setup. Internally, pyMBE uses Pint library to deal with unit transformations, which in turn should be  used by the user to define its own variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now create an instance of the ESPResSo system where we will place our molecules (a square simulation box with length = `box_l`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Box_L = 7.5*pmb.units.nm\n",
    "\n",
    "espresso_system = espressomd.System(box_l = [Box_L.to('reduced_length').magnitude]*3)\n",
    "\n",
    "print('The side of the simulation box is ', Box_L, '=' ,Box_L.to('reduced_length'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to create particles <a class=\"anchor\" id=\"particles\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particles are the smaller objects in the simulation box, which can represent small ions or other small chemical species. In turn, particles can also be used as building blocks for larger molecules or polymers, where they can represent one monomeric unit or part of it. Particle objects are used in pyMBE as input for several of its funcionalities, including to create larger molecules and peptides. The basic properties of a particle are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cation_name = 'Na'\n",
    "pmb.define_particle(name = cation_name, \n",
    "                    z = 0, \n",
    "                    sigma = 0.35*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The properties of the particles are stored in a pandas Dataframe (df). For displaying all the information of the particles one can use the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can use pyMBE to create any number of the defined particles into the ESPResSo system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_cations = 20\n",
    "pmb.create_pmb_object(name = cation_name,\n",
    "                      number_of_objects = N_cations,\n",
    "                      espresso_system = espresso_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the new set of particles..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see the particles that we have created by visualizing our ESPResSo system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picture_name = 'cation_system.png'\n",
    "do_snapshot_espresso_system(espresso_system = espresso_system, \n",
    "                               filename = picture_name)\n",
    "img = Image.open(picture_name)\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To delete an pyMBE object from the system we can use the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.destroy_pmb_object_in_system(name = cation_name, \n",
    "                                 espresso_system = espresso_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the df should be empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## How to create simple polymers <a class=\"anchor\" id=\"simple_polymers\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pyMBE can be used to easily construct coarse-grained models of simple polymers. Let us consider a coarse grained model for polydehydroalanaline (PDha) (figure below) in which its monomeric unit can be represented by three beads, as depicted in the schematics below: a backbone bead (grey), a bead for the carboxylic acid group (red) and a bead for the amino group (blue).\n",
    "\n",
    "<img src=\"../figs/PDha.png\" width=150 height=150 />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set up such polymer with pyMBE first one has to define the different particles in the monomer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDha_backbone_bead = 'BB-PDha'\n",
    "PDha_carboxyl_bead = 'COOH-PDha'\n",
    "PDha_amine_bead = 'NH3-PDha'\n",
    "\n",
    "pmb.define_particle(name = PDha_backbone_bead, \n",
    "                    z = 0, \n",
    "                    sigma = 0.4*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDha_carboxyl_bead, \n",
    "                    z = 0, \n",
    "                    sigma = 0.5*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDha_amine_bead, \n",
    "                    z = 0, \n",
    "                    sigma = 0.3*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, one defines the structure of the residue of the polymer. A residue is composed by a `central_bead` where one or various `side_chains` are attached. Each side chain can contain one particle or other residues. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDha_residue = 'PDha_mon'\n",
    "\n",
    "pmb.define_residue(name = PDha_residue, \n",
    "                   central_bead = PDha_backbone_bead,\n",
    "                   side_chains = [PDha_carboxyl_bead, PDha_amine_bead])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once done, one has to define a bond for each different type of bond in the polymer. For simplicity, in this tutorial we assume that all bonds are equal and we set-up all bonds using a harmonic potential with the following arbitrary parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bond_type = 'harmonic'\n",
    "generic_bond_lenght=0.4 * pmb.units.nm\n",
    "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n",
    "\n",
    "harmonic_bond = {'r_0'    : generic_bond_lenght,\n",
    "                 'k'      : generic_harmonic_constant\n",
    "                }\n",
    "\n",
    "pmb.define_bond(bond_type = bond_type,\n",
    "                bond_parameters = harmonic_bond,\n",
    "                particle_pairs = [[PDha_backbone_bead, PDha_backbone_bead],\n",
    "                                  [PDha_backbone_bead, PDha_carboxyl_bead],\n",
    "                                  [PDha_backbone_bead, PDha_amine_bead]])\n",
    "\n",
    "pmb.add_bonds_to_espresso(espresso_system = espresso_system)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: Currently, only harmonic and FENE bonds are supported."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, one can use the residues to define the polymer sequence given by the argument `residue_list`. One needs to add one residue in `residue_list` per each residue in the polymer chain. For instance a decamer should be created as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDha_polymer = 'PDha'\n",
    "N_monomers = 10\n",
    "\n",
    "pmb.define_molecule(name = PDha_polymer,\n",
    "                    residue_list = [PDha_residue]*N_monomers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining the polymer, we are ready to create one PdHa polymer in the center of the simulation box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_polymers = 1\n",
    "\n",
    "pmb.create_pmb_object(name = PDha_polymer, \n",
    "                      number_of_objects = N_polymers,\n",
    "                      espresso_system = espresso_system, \n",
    "                      position = [[Box_L.to('reduced_length').magnitude/2]*3]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can always track our particles..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us see what we have created..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picture_name = 'PDha_system.png'\n",
    "do_snapshot_espresso_system(espresso_system = espresso_system, \n",
    "                              filename = picture_name)\n",
    "img = Image.open(picture_name)\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete the particles and check that our df is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.destroy_pmb_object_in_system(name = PDha_polymer, \n",
    "                                 espresso_system = espresso_system)\n",
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to create complex polymers <a class=\"anchor\" id=\"complex_polymers\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pyMBE can also be used to setup models that requiere more complex side chains, i.e. with more than one bead per side chain. One example of these complex molecules is the poly(N,N-diallylglutamate) (PDAGA), whose structure is depicted in the figure below. Following the logic of the previous example, one would construct PDAGA with pyMBE by defining a `residue` with a `central_bead` for the polymer backbone (grey) and a `side_chain` attached to it. In this case, the group in the side chain of the PDAGA monomer has a complex structure. This group can be coarse-grained by defining another `residue`  composed by a new `central_bead` which represents the cyclic amine group (blue) and two `side_chains` ($\\alpha$ and $\\beta$ carboxyl) attached to it (red and orange).\n",
    "\n",
    "<img src=\"../figs/PDAGA.png\" width=150 height=150 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can start by defining each different bead of the PDAGA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDAGA_backbone_bead = 'BB-PDAGA'\n",
    "PDAGA_cyclic_amine_bead = 'NH3-PDAGA'\n",
    "PDAGA_alpha_carboxyl_bead = 'aCOOH-PDAGA'\n",
    "PDAGA_beta_carboxyl_bead = 'bCOOH-PDAGA'\n",
    "\n",
    "pmb.define_particle(name = PDAGA_backbone_bead, \n",
    "                    z = 0,\n",
    "                    sigma = 0.4*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDAGA_cyclic_amine_bead, \n",
    "                    z = 0, \n",
    "                    sigma = 0.3*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDAGA_alpha_carboxyl_bead, \n",
    "                    z = 0, \n",
    "                    sigma = 0.2*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDAGA_beta_carboxyl_bead, \n",
    "                    z = 0, \n",
    "                    sigma = 0.4*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to define the two different residues: \n",
    "1. The side chain: two carboxyl beads attached to the cyclic amine bead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDAGA_side_chain_residue = 'PDAGA_side_chain_residue'\n",
    "\n",
    "pmb.define_residue (name = PDAGA_side_chain_residue,\n",
    "                    central_bead = PDAGA_cyclic_amine_bead,\n",
    "                    side_chains = [PDAGA_alpha_carboxyl_bead, PDAGA_beta_carboxyl_bead])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Each monomeric unit of the PDAGA: the side chain defined above attached to the backbone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDAGA_monomer_residue = 'PDAGA_monomer_residue'\n",
    "pmb.define_residue( name = PDAGA_monomer_residue,\n",
    "                    central_bead = PDAGA_backbone_bead,\n",
    "                    side_chains = [PDAGA_side_chain_residue])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we need to set the bonds between the particles in a similar way as for the case of the simple polymer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bond_type = 'harmonic'\n",
    "generic_bond_lenght=0.4 * pmb.units.nm\n",
    "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n",
    "\n",
    "harmonic_bond = {'r_0'    : generic_bond_lenght,\n",
    "                 'k'      : generic_harmonic_constant,\n",
    "                 }\n",
    "\n",
    "pmb.define_bond(bond_type = bond_type,\n",
    "                bond_parameters = harmonic_bond,\n",
    "                particle_pairs = [[PDAGA_backbone_bead, PDAGA_backbone_bead],\n",
    "                                  [PDAGA_backbone_bead, PDAGA_cyclic_amine_bead],\n",
    "                                  [PDAGA_alpha_carboxyl_bead, PDAGA_cyclic_amine_bead],\n",
    "                                  [PDAGA_beta_carboxyl_bead, PDAGA_cyclic_amine_bead]])\n",
    "\n",
    "pmb.add_bonds_to_espresso(espresso_system = espresso_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us define an octamer of PDAGA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDAGA_polymer = 'PDAGA'\n",
    "N_monomers = 8\n",
    "\n",
    "pmb.define_molecule(name = PDAGA_polymer,\n",
    "                    residue_list = [PDAGA_monomer_residue]*N_monomers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we are able to create a PDAGA polymer into the ESPResSo system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_polymers = 1\n",
    "\n",
    "pmb.create_pmb_object(name = PDAGA_polymer,\n",
    "                      number_of_objects = N_polymers,\n",
    "                      espresso_system = espresso_system,\n",
    "                      position = [[Box_L.to('reduced_length').magnitude/2]*3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us see our PDAGA molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picture_name = 'PDAGA_system.png'\n",
    "do_snapshot_espresso_system(espresso_system = espresso_system, \n",
    "                               filename = picture_name)\n",
    "img = Image.open(picture_name)\n",
    "img.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete the particles and check that our df is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.destroy_pmb_object_in_system(name = PDAGA_polymer, \n",
    "                                 espresso_system = espresso_system)\n",
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to create di-block copolymers <a class=\"anchor\" id=\"diblock_copolymers\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In turn, the residues previously defined to build the PDAGA and PDha molecules can be used to build more complex polymers such as a di-block PDha-PDAGA copolymer, as shown in the picture below\n",
    "\n",
    "<img src=\"../figs/PDAGA_PDha_diblock_copolymer.png\" width=250 height=250 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining each different bead of the PDha and PDAGA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PDha_backbone_bead = 'BB-PDha'\n",
    "PDha_carboxyl_bead = 'COOH-PDha'\n",
    "PDha_amine_bead = 'NH3-PDha'\n",
    "\n",
    "pmb.define_particle(name = PDha_backbone_bead, \n",
    "                    z = 0,\n",
    "                    sigma = 0.4*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDha_carboxyl_bead,\n",
    "                    z = 0,\n",
    "                    sigma = 0.5*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDha_amine_bead,\n",
    "                    z = 0,\n",
    "                    sigma = 0.3*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "PDAGA_backbone_bead = 'BB-PDAGA'\n",
    "PDAGA_cyclic_amine_bead = 'NH3-PDAGA'\n",
    "PDAGA_alpha_carboxyl_bead = 'aCOOH-PDAGA'\n",
    "PDAGA_beta_carboxyl_bead = 'bCOOH-PDAGA'\n",
    "\n",
    "\n",
    "pmb.define_particle(name = PDAGA_backbone_bead,\n",
    "                    z = 0,\n",
    "                    sigma = 0.4*pmb.units.nm, \n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDAGA_cyclic_amine_bead,\n",
    "                    z = 0, \n",
    "                    sigma = 0.3*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDAGA_alpha_carboxyl_bead,\n",
    "                    z = 0,\n",
    "                    sigma = 0.2*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))\n",
    "\n",
    "pmb.define_particle(name = PDAGA_beta_carboxyl_bead, \n",
    "                    z = 0,\n",
    "                    sigma = 0.4*pmb.units.nm,\n",
    "                    epsilon = 1*pmb.units('reduced_energy'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining the different residues for the PDha and PDAGA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bond_type = 'harmonic'\n",
    "generic_bond_lenght=0.4 * pmb.units.nm\n",
    "generic_harmonic_constant = 400 * pmb.units('reduced_energy / reduced_length**2')\n",
    "\n",
    "harmonic_bond = {'r_0'    : generic_bond_lenght,\n",
    "                 'k'      : generic_harmonic_constant,\n",
    "                 }\n",
    "\n",
    "######################################################################\n",
    "################################# PDha ###############################\n",
    "######################################################################\n",
    "\n",
    "PDha_residue = 'PDha_residue'\n",
    "pmb.define_residue(name = PDha_residue, \n",
    "                   central_bead =  PDha_backbone_bead ,\n",
    "                   side_chains = [PDha_carboxyl_bead,PDha_amine_bead])\n",
    "\n",
    "\n",
    "pmb.define_bond(bond_type = bond_type,\n",
    "                bond_parameters = harmonic_bond,\n",
    "                particle_pairs = [[PDha_backbone_bead, PDha_backbone_bead],\n",
    "                                  [PDha_backbone_bead, PDha_carboxyl_bead],\n",
    "                                  [PDha_backbone_bead, PDha_amine_bead]])\n",
    "\n",
    "######################################################################\n",
    "################################ PDAGA ###############################\n",
    "######################################################################\n",
    "\n",
    "PDAGA_monomer_residue = 'PDAGA_monomer_residue'\n",
    "pmb.define_residue( name = PDAGA_monomer_residue,\n",
    "                    central_bead = PDAGA_backbone_bead,\n",
    "                    side_chains = [PDAGA_side_chain_residue])\n",
    "\n",
    "PDAGA_side_chain_residue = 'PDAGA_side_chain_residue'\n",
    "pmb.define_residue (name = PDAGA_side_chain_residue,\n",
    "                    central_bead = PDAGA_cyclic_amine_bead,\n",
    "                    side_chains = [PDAGA_alpha_carboxyl_bead,PDAGA_beta_carboxyl_bead])\n",
    "\n",
    "pmb.define_bond(bond_type = bond_type,\n",
    "                bond_parameters = harmonic_bond,\n",
    "                particle_pairs = [[PDAGA_backbone_bead, PDAGA_backbone_bead],\n",
    "                                  [PDAGA_backbone_bead, PDAGA_cyclic_amine_bead],\n",
    "                                  [PDAGA_alpha_carboxyl_bead, PDAGA_cyclic_amine_bead],\n",
    "                                  [PDAGA_beta_carboxyl_bead, PDAGA_cyclic_amine_bead]])\n",
    "\n",
    "######################################################################\n",
    "############################# PDha - PDAGA ###########################\n",
    "######################################################################\n",
    "\n",
    "pmb.define_bond(bond_type = bond_type,\n",
    "                bond_parameters = harmonic_bond,\n",
    "                particle_pairs = [[PDha_backbone_bead, PDAGA_backbone_bead]])\n",
    "\n",
    "pmb.add_bonds_to_espresso(espresso_system = espresso_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining the di-block polymer molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_monomers_PDha = 4\n",
    "N_monomers_PDAGA = 4\n",
    "diblock_polymer = 'diblock'\n",
    "\n",
    "pmb.define_molecule(name = diblock_polymer,\n",
    "                    residue_list = [PDha_residue]*N_monomers_PDha+[PDAGA_monomer_residue]*N_monomers_PDAGA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating the di-block polymer into the ESPResSo system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_polymers = 1\n",
    "\n",
    "pmb.create_pmb_object(name = diblock_polymer,\n",
    "                      number_of_objects = N_polymers,\n",
    "                      espresso_system = espresso_system,\n",
    "                      position = [[Box_L.to('reduced_length').magnitude/2]*3]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us see our di-block PDha-PDAGA molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picture_name = 'diblock_system.png'\n",
    "do_snapshot_espresso_system(espresso_system = espresso_system, \n",
    "                               filename = picture_name)\n",
    "img = Image.open(picture_name)\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete the particles and check that our df is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.destroy_pmb_object_in_system(name = diblock_polymer, \n",
    "                     espresso_system = espresso_system)\n",
    "pmb.filter_df(pmb_type = 'particle')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Practice by creating an alternating copolymer <a class=\"anchor\" id=\"alternating_copolymers\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, one can set-up an alternating PAA-PVAm copolymer (Polyacrylic acid - Polyvinylamine), whose structure is depicted in the figure below.\n",
    "\n",
    "<img src=\"../figs/PAA_PVAm.png\" width=250 height=250 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create an alternating copolymer one needs to set up the system in a similar way as we did it for the di-block copolymer but you should define differently the 'residue_list' in the function 'pmb.define_molecule', as follows:\n",
    "\n",
    "residue_list = [ 'Residue_1' , 'Residue_2' ] * 'N_monomers'\n",
    "\n",
    "Let's practice by creating the alternating copolymer (PAA-PVAm)$_{4}$.\n",
    "\n",
    "### Tasks to do:\n",
    "\n",
    "1. Define each different bead of the PAA and PVAm.\n",
    "2. Define the different residues for the PAA and PVAm.\n",
    "3. Define the alternating block copolymer molecule. \n",
    "4. Create the alternating block copolymer into the ESPResSo system.\n",
    "5. Visualize your creation.\n",
    "6. Delete the molecule and check that your df is empty.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Define each different bead of the PAA and PVAm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Define the different residues for the PAA and PVAm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Define the alternating block copolymer molecule. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. Create the alternating block copolymer into the ESPResSo system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5. Visualize your creation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 6. Delete the molecule and check that your df is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to create peptides <a class=\"anchor\" id=\"peptides\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pyMBE includes built-on functions to facilitate the setting up of coarse-grained models for peptides from their aminoacid sequence. Currently, there are two different coarse-grained models implemented: \n",
    "\n",
    "* `1beadAA`, where the aminoacid is represented by one single bead.\n",
    "* `2beadAA`, where the aminoacid is represented by two beads (backbone and side-chain). \n",
    "\n",
    "We provide reference parameters in the folder (`parameters`) which can be loaded into pyMBE. The peptide sequence should be provided as a str composed either by the list of the one letter code or the list of the three letter code of the corresponding aminoacids. For example, the two possible ways to provide the peptide Cysteine$_3$ - Glutamic acid$_2$ - Histidine$_4$ - Valine are:\n",
    "\n",
    "* one letter code: 'CCCEEHHHHV'\n",
    "* three letter code: 'CYS-CYS-CYS-GLU-GLU-HIS-HIS-HIS-HIS-VAL'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set up the peptide Lysine$_5$ - Glutamic acid$_5$ using a two beads coarse-grained model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_peptide = 1\n",
    "sequence = \"KKKKKEEEEE\"\n",
    "model = '2beadAA'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the peptide parametrization reported by Lunkad et al., which is provided in the reference folder. After loading the parameters we should add the bonds to the ESPResSo system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.load_interaction_parameters(filename = pmb.get_resource('parameters/peptides/Lunkad2021.json'))\n",
    "pmb.add_bonds_to_espresso(espresso_system = espresso_system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pmb.df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can define our peptide and create it into the ESPResSo system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.define_peptide(name = sequence, \n",
    "                   sequence = sequence, \n",
    "                   model = model)\n",
    "\n",
    "pmb.create_pmb_object(name = sequence,\n",
    "                      number_of_objects = N_peptide,\n",
    "                      espresso_system = espresso_system,\n",
    "                      position = [[Box_L.to('reduced_length').magnitude/2]*3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us visualize our peptide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picture_name = 'alternating_system.png'\n",
    "do_snapshot_espresso_system(espresso_system = espresso_system, \n",
    "                               filename = picture_name)\n",
    "img = Image.open(picture_name)\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Delete the particles and check that our df is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmb.destroy_pmb_object_in_system(name = sequence, \n",
    "                                 espresso_system = espresso_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the setup using the three letter code. pyMBE automatically detects and transforms into the one-letter code using its own protein sequence parser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequence = 'LYS-LYS-LYS-LYS-LYS-GLU-GLU-GLU-GLU-GLU'\n",
    "\n",
    "pmb.define_peptide(name = sequence, \n",
    "                   sequence = sequence, \n",
    "                   model = model)\n",
    "\n",
    "print('one letter code', pmb.protein_sequence_parser(sequence=sequence))\n",
    "print('defined peptide sequence ', sequence)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lunkad, R. et al.  Molecular Systems Design & Engineering (2021), 6(2), 122-131."
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "pyMBE",
   "language": "python",
   "name": "pymbe"
  },
  "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.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}