{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 07. Molecules CG Mapping \n", "\n", "In this tutorial, we show how to generate a CG mapping matrix for a molecule given a bead mapping. The trajectory and topology file come from a AA simulation done in gromacs (see `Molecules_CG_Mapping` folder). The protein is FF (diphenylalanine) and the solvent is a mixture of water and methanol." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# disable GPU\n", "import os\n", "os.environ['CUDA_VISIBLE_DEVICES'] = '-1'\n", "import hoomd, hoomd.htf as htf, hoomd.md\n", "import matplotlib.pyplot as plt\n", "import tensorflow as tf\n", "import numpy as np\n", "import MDAnalysis as mda\n", "from os import path\n", "import matplotlib.pyplot as plt,matplotlib" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/mgholiza/.conda/envs/hoomd-tf2/lib/python3.7/site-packages/MDAnalysis/core/universe.py:171: UserWarning: No coordinate reader found for Molecules_CG_Mapping/nvt_prod.tpr. Skipping this file.\n", " 'this file.'.format(filename))\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "[[0.8224383 0.05918723 0.05918723 0.05918723 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0.44409524 0.03726984\n", " 0.44409524 0.03726984 0.03726984 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0.15577257 0.15577257 0.01307291\n", " 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291\n", " 0.15577257 0.01307291 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0.42880501 0.57119499 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0.93286579 0.06713421\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0.44409524 0.03726984 0.44409524 0.03726984 0.03726984 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.15577257\n", " 0.15577257 0.01307291 0.15577257 0.01307291 0.15577257 0.01307291\n", " 0.15577257 0.01307291 0.15577257 0.01307291 0. 0.\n", " 0. ]\n", " [0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0.27291648 0.36354176\n", " 0.36354176]]\n" ] } ], "source": [ "# Loading inputs\n", "TPR = 'Molecules_CG_Mapping/nvt_prod.tpr'\n", "tpr = mda.Universe(TPR)\n", "TRAJECTORY = 'Molecules_CG_Mapping/traj.trr'\n", "u = mda.Universe(TPR, TRAJECTORY)\n", "\n", "# Generating Mapping Matrix for FF\n", "protein_FF = u.select_atoms(\"resname PHE and resid 0:1\")\n", "Beads_distribution = [['N','H1','H2','H3'],\n", " ['CA','HA','CB','HB1','HB2'],\n", " ['CG','CD1','HD1','CD2','HD2','CE1','HE1','CE2','HE2','CZ','HZ'],\n", " ['C','O'],\n", " ['N','H'],\n", " ['CA','HA','CB','HB1','HB2'],\n", " ['CG','CD1','HD1','CD2','HD2','CE1','HE1','CE2','HE2','CZ','HZ'],\n", " ['C','O1','O2']]\n", "mapping_FF = htf.matrix_mapping(protein_FF,Beads_distribution)\n", "print (mapping_FF)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.88809574 0.05595213 0.05595213]]\n" ] } ], "source": [ "# Generating Mapping Matrix for Water\n", "water = u.select_atoms(\"resname SOL and resid 500\")\n", "Beads_distribution = [['OW','HW1','HW2']]\n", "mapping_water = htf.matrix_mapping(water,Beads_distribution) \n", "print (mapping_water)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.37484707 0.03145832 0.03145832 0.03145832 0.49931966 0.03145832]]\n" ] } ], "source": [ "# Generating Mapping Matrix for Methanol\n", "methanol = u.select_atoms(\"resname MET and resid 11665 \")\n", "Beads_distribution_methanol = [['C','H','H','H','OA','HO']]\n", "mapping_methanol = htf.matrix_mapping(methanol,Beads_distribution_methanol) \n", "print (mapping_methanol)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Getting the segment id of each molecule in topology\n", "_,idx = np.unique(u.select_atoms('all').segids,return_index=True)\n", "seg_id_list = u.select_atoms('all').segids[np.sort(idx)].tolist()\n", "\n", "# Getting the list of every molecule type name in topology\n", "_,idx = np.unique(u.atoms.resnames,return_index=True)\n", "resname_list = u.atoms.resnames[np.sort(idx)].tolist()\n", "\n", "# Getting list of atoms in each type of molecule\n", "atoms_in_molecule_list = [protein_FF.names,\n", " water.names,\n", " methanol.names]" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "160\n" ] } ], "source": [ "atoms_in_molecule_list = [protein_FF.names]\n", "molecule_mapping_index = htf.find_molecules_from_topology(u,atoms_in_molecule_list, selection = \"resname PHE\")\n", "molecule_mapping = mapping_FF\n", "\n", "# get number of atoms\n", "N = sum([len(m) for m in molecule_mapping_index])\n", "# get number of molecules\n", "M = len(molecule_mapping_index)\n", "# get number of beads\n", "B = molecule_mapping.shape[0]\n", "\n", "#create a mass-weighted (M * bead_number) x N mapping operator \n", "cg_mapping = htf.sparse_mapping([molecule_mapping for _ in molecule_mapping_index], \n", " molecule_mapping_index)\n", "assert cg_mapping.shape == (M * B, N)\n", "print (M*B)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }