Skip to content

Instantly share code, notes, and snippets.

@yuxuanzhuang
Last active February 2, 2021 14:29
Show Gist options
  • Save yuxuanzhuang/d765797bc876b8e731549e82ba7547f5 to your computer and use it in GitHub Desktop.
Save yuxuanzhuang/d765797bc876b8e731549e82ba7547f5 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"# Modified from penetest.py file in CHARMM-GUI"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import argparse\n",
"import sys\n",
"import networkx as nx\n",
"import re\n",
"import numpy as np\n",
"import itertools"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import MDAnalysis as mda"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"_s = re.compile('\\s+')\n",
"_p = re.compile('(\\d+)\\s+(\\d+)')"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def lsqp(atoms):\n",
" com = atoms.mean(axis=0)\n",
" #u, d, v = np.linalg.svd(atoms-com)\n",
"\n",
" axes = np.zeros((len(atoms), 3))\n",
" for i in range(len(atoms)):\n",
" p1 = atoms[i]\n",
" if i == len(atoms)-1:\n",
" p2 = atoms[0]\n",
" else:\n",
" p2 = atoms[i+1]\n",
" a = np.cross(p1, p2)\n",
" axes += a\n",
" u, d, v = np.linalg.svd(axes)\n",
" i = 0\n",
" d = -np.dot(v[i], com)\n",
" n = -np.array((v[i,0], v[i,1], d))/v[i,2]\n",
" return v[i], com, n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def intriangle(triangle, axis, u, p):\n",
" # http://www.softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm\n",
" p1, p2, p3 = triangle\n",
" w0 = p - p1\n",
" a = -np.dot(axis, w0)\n",
" b = np.dot(axis, u)\n",
" if (abs(b) < 0.01): return False\n",
"\n",
" r = a / b\n",
" if r < 0.0: return False\n",
" if r > 1.0: return False\n",
"\n",
" I = p + u * r\n",
"\n",
" u = p2 - p1\n",
" v = p3 - p1\n",
" uu = np.dot(u, u)\n",
" uv = np.dot(u, v)\n",
" vv = np.dot(v, v)\n",
" w = I - p1\n",
" wu = np.dot(w, u)\n",
" wv = np.dot(w, v)\n",
" D = uv * uv - uu * vv\n",
"\n",
" s = (uv * wv - vv * wu)/D\n",
" if (s < 0 or s > 1): return False\n",
" t = (uv * wu - uu * wv)/D\n",
" if (t < 0 or (s+t) > 1): return False\n",
" return True"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"def build_topology(universe, selection):\n",
" g = nx.Graph()\n",
"\n",
" # Atoms\n",
" natom = universe.atoms.n_atoms\n",
" for atom in universe.select_atoms(selection).atoms: # might be buggy\n",
" g.add_node(atom.index + 1, **{'segid': atom.segid,\n",
" 'resname': atom.resname,\n",
" 'name': atom.name,\n",
" 'resid': atom.resid})\n",
" # Bonds\n",
" for bond in universe.select_atoms(selection).bonds:\n",
" num1, num2 = bond.atoms.indices + 1\n",
" if g.has_node(num1) and g.has_node(num2):\n",
" g.add_edge(num1, num2)\n",
" return g"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"def check_ring_penetration(top, coord, pbc=[], xtl='rect', verbose=0):\n",
" # ring penetration test\n",
" # 1. find rings\n",
" # 2. build least square plane\n",
" # 3. project atoms ring constituent atoms onto the plane and build convex\n",
" # 4. find two bonded atoms that are at the opposite side of the plane\n",
" # 5. determine the point of intersection is enclosed in the ring\n",
" #\n",
" from networkx.algorithms.components import connected_components\n",
" molecules = (top.subgraph(c) for c in connected_components(top))\n",
"\n",
" allatoms = np.array([coord[num] for num in top.nodes()])\n",
" atoms_map = np.array([num for num in top.nodes()])\n",
" natoms = len(allatoms)\n",
" if pbc:\n",
" atoms_map_reverse = {}\n",
" for i,num in enumerate(top.nodes()):\n",
" atoms_map_reverse[num] = i\n",
"\n",
" a = float(pbc[0])\n",
" b = float(pbc[1])\n",
" n = len(allatoms)\n",
" if xtl == 'rect':\n",
" allatoms = np.tile(allatoms, (9,1))\n",
" op = ((a,0),(a,b),(0,b),(-a,b),(-a,0),(-a,-b),(0,-b),(a,-b))\n",
" for i in range(8):\n",
" x,y = op[i]\n",
" allatoms[n*(i+1):n*(i+2),0] += x\n",
" allatoms[n*(i+1):n*(i+2),1] += y\n",
" atoms_map = np.tile(atoms_map, 9)\n",
" if xtl =='hexa':\n",
" allatoms = np.tile(allatoms, (7,1))\n",
" rot = lambda theta: np.matrix(((np.cos(np.radians(theta)), -np.sin(np.radians(theta))),\n",
" (np.sin(np.radians(theta)), np.cos(np.radians(theta)))))\n",
" op = (rot(15), rot(75), rot(135), rot(195), rot(255), rot(315))\n",
" d = np.array((a, 0))\n",
" for i in range(6):\n",
" xy = np.dot(d, op[i])\n",
" allatoms[n*(i+1):n*(i+2),:2] = allatoms[n*(i+1):n*(i+2),:2] + xy\n",
" atoms_map = np.tile(atoms_map, 7)\n",
"\n",
" # print out image atoms\n",
" #fp = open('image.pdb', 'w')\n",
" #for i,atom in enumerate(allatoms):\n",
" # x, y, z = atom\n",
" # fp.write(\"HETATM%5d %-3s %3s %4d %8.3f%8.3f%8.3f 0.00 0.00 \\n\" % (i, 'C', 'DUM', i, x, y, z))\n",
"\n",
" pen_pairs = []\n",
" pen_cycles = []\n",
"\n",
" for m in molecules:\n",
" cycles = nx.cycle_basis(m)\n",
" if not cycles: continue\n",
" for cycle in cycles:\n",
" flag = False\n",
" atoms = np.array([coord[num] for num in cycle])\n",
" if len(set([top.nodes[num]['resid'] for num in cycle])) > 1: continue\n",
" if verbose:\n",
" num = cycle[0]\n",
" print('found ring:', top.nodes[num]['segid'], top.nodes[num]['resid'], top.nodes[num]['resname'])\n",
"\n",
" # build least square fit plane\n",
" axis, com, n = lsqp(atoms)\n",
"\n",
" # project atoms to the least square fit plane\n",
" for i,atom in enumerate(atoms):\n",
" w = np.dot(axis, atom-com)*axis + com\n",
" atoms[i] = com + (atom - w)\n",
"\n",
" maxd = np.max(np.sqrt(np.sum(np.square(atoms - com), axis=1)))\n",
"\n",
" d = np.sqrt(np.sum(np.square(allatoms-com), axis=1))\n",
" nums = np.squeeze(np.argwhere(d < 3))\n",
"\n",
" # find two bonded atoms that are at the opposite side of the plane\n",
" for num in nums:\n",
" num1 = atoms_map[num]\n",
"\n",
" for num2 in top[num1]:\n",
" if num1 in cycle or num2 in cycle: continue\n",
" if num > natoms:\n",
" # image atoms\n",
" offset = int(num / natoms)\n",
" coord1 = allatoms[num]\n",
" coord2 = allatoms[atoms_map_reverse[num2] + offset * natoms]\n",
" else:\n",
" coord1 = coord[num1]\n",
" coord2 = coord[num2]\n",
"\n",
" v1 = np.dot(coord1 - com, axis)\n",
" v2 = np.dot(coord2 - com, axis)\n",
" if v1 * v2 > 0: continue\n",
"\n",
" # point of intersection of the least square fit plane\n",
" s = -np.dot(axis, coord1-com)/np.dot(axis, coord2-coord1)\n",
" p = coord1 + s*(coord2-coord1)\n",
"\n",
" d = np.sqrt(np.sum(np.square(p-com)))\n",
" if d > maxd: continue\n",
" if verbose:\n",
" print('found potentially pentrarting bond:',\n",
" top.nodes[num1]['segid'],\n",
" top.nodes[num1]['resid'],\n",
" top.nodes[num1]['resname'],\n",
" top.nodes[num1]['name'],\n",
" top.nodes[num2]['name'])\n",
"\n",
" d = 0\n",
" for i in range(0, len(atoms)):\n",
" p1 = atoms[i] - p\n",
" try: p2 = atoms[i+1] - p\n",
" except: p2 = atoms[0] - p\n",
" d += np.arccos(np.dot(p1, p2)/np.linalg.norm(p1)/np.linalg.norm(p2))\n",
"\n",
" wn = d/2/np.pi\n",
" if wn > 0.9 and wn < 1.1:\n",
" # we have a case\n",
" pen_pairs.append((num1, num2))\n",
" pen_cycles.append(cycle)\n",
" flag = True\n",
" break\n",
"\n",
" if flag: break\n",
"\n",
" return pen_pairs, pen_cycles"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"def check_universe_ring_penetration(universe):\n",
" selection = 'not resname TIP3 and not (name H*)'\n",
" \n",
" top = build_topology(universe, selection)\n",
" ag = universe.select_atoms(selection)\n",
" for frame, ts in enumerate(universe.trajectory):\n",
" print('In frame', frame)\n",
" coord = dict(zip(ag.ids + 1, ag.positions))\n",
" if len(top.nodes()) != len(coord):\n",
" raise AtomMismatch('Number of atoms does not match')\n",
"\n",
" print('%lipid ring pentration')\n",
" \n",
" # only rect pbc have been tested\n",
" pairs, rings = check_ring_penetration(top, coord)\n",
" if pairs:\n",
" print('found a ring penetration:')\n",
" for i, cycle in enumerate(rings):\n",
" print('- %s %s %s %s | %s %s %s %s' % (top.nodes[pairs[i][0]]['segid'], top.nodes[pairs[i][0]]['resid'], top.nodes[pairs[i][0]]['resname'], ' '.join([top.nodes[num]['name'] for num in pairs[i]]), top.nodes[cycle[0]]['segid'], top.nodes[cycle[0]]['resid'], top.nodes[cycle[0]]['resname'], ' '.join([top.nodes[num]['name'] for num in cycle])))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In frame 0\n",
"%lipid ring pentration\n",
"found a ring penetration:\n",
"- seg_7_SAPI24 2779 SAPI24 C36 C37 | seg_0_PROA 946 TRP CE2 CZ2 CH2 CZ3 CE3 CD2\n",
"- seg_5_PLPE 2398 PLPE C23 C24 | seg_0_PROA 1266 TRP CD1 NE1 CE2 CD2 CG\n",
"- seg_3_PLPI 2123 PLPI C34 C35 | seg_0_PROA 1264 PRO CD CG CB CA N\n",
"- seg_7_SAPI24 2772 SAPI24 C38 C39 | seg_1_PROB 1665 PHE CD1 CE1 CZ CE2 CD2 CG\n",
"- seg_4_POPC 2261 POPC C23 C24 | seg_2_CHL1 1947 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_3_PLPI 2125 PLPI C310 C311 | seg_2_CHL1 1956 CHL1 C14 C13 C12 C11 C9 C8\n",
"- seg_5_PLPE 2454 PLPE C34 C35 | seg_2_CHL1 1970 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_POPC 2331 POPC C213 C214 | seg_2_CHL1 1978 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_4_POPC 2372 POPC C36 C37 | seg_2_CHL1 1982 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2623 PLPC C215 C216 | seg_2_CHL1 1983 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_9_PLPA 2881 PLPA C310 C311 | seg_2_CHL1 1992 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_9_PLPA 2879 PLPA C212 C213 | seg_2_CHL1 1997 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2604 PLPC C212 C213 | seg_2_CHL1 1998 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_3_PLPI 2166 PLPI C311 C312 | seg_2_CHL1 2001 CHL1 C6 C7 C8 C9 C10 C5\n",
"- seg_6_PLPC 2523 PLPC C26 C27 | seg_2_CHL1 2009 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_6_PLPC 2523 PLPC C215 C216 | seg_2_CHL1 2009 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2601 PLPC C28 C29 | seg_2_CHL1 2014 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_8_POPA 2825 POPA C25 C26 | seg_2_CHL1 2017 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_9_PLPA 2869 PLPA C26 C27 | seg_2_CHL1 2019 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2589 PLPC C310 C311 | seg_2_CHL1 2027 CHL1 C14 C15 C16 C17 C13\n",
"- seg_6_PLPC 2564 PLPC C315 C316 | seg_2_CHL1 2027 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2532 PLPC C216 C217 | seg_2_CHL1 2028 CHL1 C8 C7 C6 C5 C10 C9\n",
"- seg_5_PLPE 2404 PLPE C313 C314 | seg_2_CHL1 2036 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_POPC 2270 POPC O21 C21 | seg_2_CHL1 2041 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_3_PLPI 2186 PLPI C315 C316 | seg_2_CHL1 2042 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_9_PLPA 2861 PLPA C25 C26 | seg_2_CHL1 2048 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_6_PLPC 2598 PLPC C210 C211 | seg_2_CHL1 2062 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_3_PLPI 2123 PLPI C215 C216 | seg_2_CHL1 2066 CHL1 C16 C15 C14 C13 C17\n",
"- seg_6_PLPC 2564 PLPC C213 C214 | seg_2_CHL1 2069 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_4_POPC 2334 POPC C35 C36 | seg_2_CHL1 2070 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_4_POPC 2334 POPC C39 C310 | seg_2_CHL1 2070 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2637 PLPC C39 C310 | seg_2_CHL1 2084 CHL1 C6 C5 C10 C9 C8 C7\n",
"- seg_6_PLPC 2569 PLPC C36 C37 | seg_2_CHL1 2089 CHL1 C8 C7 C6 C5 C10 C9\n",
"- seg_6_PLPC 2541 PLPC C211 C212 | seg_2_CHL1 2112 CHL1 C5 C6 C7 C8 C9 C10\n"
]
}
],
"source": [
"universe = mda.Universe('./SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/em.tpr',\n",
" './SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/em.gro')\n",
"check_universe_ring_penetration(universe)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In frame 0\n",
"%lipid ring pentration\n",
"found a ring penetration:\n",
"- seg_7_SAPI24 2772 SAPI24 C38 C39 | seg_1_PROB 1665 PHE CD1 CE1 CZ CE2 CD2 CG\n",
"- seg_4_POPC 2261 POPC C23 C24 | seg_2_CHL1 1947 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_3_PLPI 2125 PLPI C310 C311 | seg_2_CHL1 1956 CHL1 C14 C13 C12 C11 C9 C8\n",
"- seg_5_PLPE 2454 PLPE C34 C35 | seg_2_CHL1 1970 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_POPC 2331 POPC C213 C214 | seg_2_CHL1 1978 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_4_POPC 2372 POPC C36 C37 | seg_2_CHL1 1982 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2623 PLPC C215 C216 | seg_2_CHL1 1983 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_9_PLPA 2881 PLPA C310 C311 | seg_2_CHL1 1992 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_9_PLPA 2879 PLPA C212 C213 | seg_2_CHL1 1997 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2604 PLPC C212 C213 | seg_2_CHL1 1998 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_3_PLPI 2166 PLPI C311 C312 | seg_2_CHL1 2001 CHL1 C6 C7 C8 C9 C10 C5\n",
"- seg_6_PLPC 2523 PLPC C26 C27 | seg_2_CHL1 2009 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_6_PLPC 2523 PLPC C215 C216 | seg_2_CHL1 2009 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2601 PLPC C28 C29 | seg_2_CHL1 2014 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_8_POPA 2825 POPA C25 C26 | seg_2_CHL1 2017 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_9_PLPA 2869 PLPA C26 C27 | seg_2_CHL1 2019 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2589 PLPC C310 C311 | seg_2_CHL1 2027 CHL1 C14 C15 C16 C17 C13\n",
"- seg_6_PLPC 2564 PLPC C315 C316 | seg_2_CHL1 2027 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2532 PLPC C216 C217 | seg_2_CHL1 2028 CHL1 C8 C7 C6 C5 C10 C9\n",
"- seg_5_PLPE 2404 PLPE C313 C314 | seg_2_CHL1 2036 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_POPC 2270 POPC O21 C21 | seg_2_CHL1 2041 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_3_PLPI 2186 PLPI C315 C316 | seg_2_CHL1 2042 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_9_PLPA 2861 PLPA C25 C26 | seg_2_CHL1 2048 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_6_PLPC 2598 PLPC C210 C211 | seg_2_CHL1 2062 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_3_PLPI 2123 PLPI C215 C216 | seg_2_CHL1 2066 CHL1 C16 C15 C14 C13 C17\n",
"- seg_6_PLPC 2564 PLPC C213 C214 | seg_2_CHL1 2069 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_4_POPC 2334 POPC C35 C36 | seg_2_CHL1 2070 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_4_POPC 2334 POPC C39 C310 | seg_2_CHL1 2070 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2637 PLPC C39 C310 | seg_2_CHL1 2084 CHL1 C6 C5 C10 C9 C8 C7\n",
"- seg_6_PLPC 2569 PLPC C36 C37 | seg_2_CHL1 2089 CHL1 C8 C7 C6 C5 C10 C9\n",
"- seg_6_PLPC 2541 PLPC C211 C212 | seg_2_CHL1 2112 CHL1 C5 C6 C7 C8 C9 C10\n"
]
}
],
"source": [
"universe = mda.Universe('./SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/alchem.tpr',\n",
" './SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/alchem.gro')\n",
"check_universe_ring_penetration(universe)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In frame 0\n",
"%lipid ring pentration\n",
"found a ring penetration:\n",
"- seg_4_POPC 2261 POPC C23 C24 | seg_2_CHL1 1947 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_3_PLPI 2125 PLPI C310 C311 | seg_2_CHL1 1956 CHL1 C14 C13 C12 C11 C9 C8\n",
"- seg_5_PLPE 2454 PLPE C34 C35 | seg_2_CHL1 1970 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_POPC 2331 POPC C213 C214 | seg_2_CHL1 1978 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_4_POPC 2372 POPC C36 C37 | seg_2_CHL1 1982 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2623 PLPC C215 C216 | seg_2_CHL1 1983 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_9_PLPA 2881 PLPA C310 C311 | seg_2_CHL1 1992 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_9_PLPA 2879 PLPA C212 C213 | seg_2_CHL1 1997 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2604 PLPC C212 C213 | seg_2_CHL1 1998 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_3_PLPI 2166 PLPI C311 C312 | seg_2_CHL1 2001 CHL1 C6 C7 C8 C9 C10 C5\n",
"- seg_6_PLPC 2523 PLPC C26 C27 | seg_2_CHL1 2009 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_6_PLPC 2523 PLPC C215 C216 | seg_2_CHL1 2009 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2601 PLPC C28 C29 | seg_2_CHL1 2014 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_8_POPA 2825 POPA C25 C26 | seg_2_CHL1 2017 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_9_PLPA 2869 PLPA C26 C27 | seg_2_CHL1 2019 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2589 PLPC C310 C311 | seg_2_CHL1 2027 CHL1 C14 C15 C16 C17 C13\n",
"- seg_6_PLPC 2564 PLPC C315 C316 | seg_2_CHL1 2027 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_6_PLPC 2532 PLPC C216 C217 | seg_2_CHL1 2028 CHL1 C8 C7 C6 C5 C10 C9\n",
"- seg_5_PLPE 2404 PLPE C313 C314 | seg_2_CHL1 2036 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_POPC 2270 POPC O21 C21 | seg_2_CHL1 2041 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_3_PLPI 2186 PLPI C315 C316 | seg_2_CHL1 2042 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_9_PLPA 2861 PLPA C25 C26 | seg_2_CHL1 2048 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_6_PLPC 2598 PLPC C210 C211 | seg_2_CHL1 2062 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_3_PLPI 2123 PLPI C215 C216 | seg_2_CHL1 2066 CHL1 C16 C15 C14 C13 C17\n",
"- seg_6_PLPC 2564 PLPC C213 C214 | seg_2_CHL1 2069 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_4_POPC 2334 POPC C35 C36 | seg_2_CHL1 2070 CHL1 C5 C4 C3 C2 C1 C10\n",
"- seg_4_POPC 2334 POPC C39 C310 | seg_2_CHL1 2070 CHL1 C14 C8 C9 C11 C12 C13\n",
"- seg_6_PLPC 2637 PLPC C39 C310 | seg_2_CHL1 2084 CHL1 C6 C5 C10 C9 C8 C7\n",
"- seg_6_PLPC 2569 PLPC C36 C37 | seg_2_CHL1 2089 CHL1 C8 C7 C6 C5 C10 C9\n",
"- seg_6_PLPC 2541 PLPC C211 C212 | seg_2_CHL1 2112 CHL1 C5 C6 C7 C8 C9 C10\n"
]
}
],
"source": [
"universe = mda.Universe('./SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/alchem_prob.tpr',\n",
" './SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/alchem_prob.gro')\n",
"check_universe_ring_penetration(universe)"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In frame 0\n",
"%lipid ring pentration\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-5-5b658d4b0d49>:17: RuntimeWarning: divide by zero encountered in true_divide\n",
" n = -np.array((v[i,0], v[i,1], d))/v[i,2]\n"
]
}
],
"source": [
"universe = mda.Universe('./SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/alchem_chl1.tpr',\n",
" './SIMULATIONS/NACHRA7_NOPNU_EPJ_SLIPIDS/alchembed/alchem_chl1.gro')\n",
"check_universe_ring_penetration(universe)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In frame 0\n",
"%lipid ring pentration\n",
"found a ring penetration:\n",
"- seg_16_PNU 153732 PNU C6 C12 | seg_4_CHL1 1970 CHL1 C14 C15 C16 C17 C13\n",
"- seg_16_PNU 153732 PNU C1 N1 | seg_4_CHL1 1970 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_16_PNU 153732 PNU C2 N2 | seg_4_CHL1 1970 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_16_PNU 153732 PNU O2 C3 | seg_4_CHL1 1970 CHL1 C3 C4 C5 C10 C1 C2\n",
"- seg_16_PNU 153734 PNU C1 N1 | seg_4_CHL1 1980 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_16_PNU 153728 PNU C1 N1 | seg_4_CHL1 2042 CHL1 C8 C14 C13 C12 C11 C9\n",
"- seg_16_PNU 153728 PNU C2 C5 | seg_4_CHL1 2042 CHL1 C5 C6 C7 C8 C9 C10\n",
"- seg_16_PNU 153728 PNU N2 O2 | seg_4_CHL1 2042 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_16_PNU 153734 PNU C3 C4 | seg_4_CHL1 2123 CHL1 C4 C5 C10 C1 C2 C3\n",
"- seg_4_CHL1 2042 CHL1 C5 C10 | seg_16_PNU 153728 PNU N2 O2 C3 C5 C2\n",
"- seg_4_CHL1 1970 CHL1 C15 C16 | seg_16_PNU 153732 PNU C7 C8 C9 C11 C12 C6\n",
"- seg_4_CHL1 1970 CHL1 C5 C10 | seg_16_PNU 153732 PNU N2 O2 C3 C5 C2\n"
]
}
],
"source": [
"universe = mda.Universe('./SIMULATIONS/NACHRA7_PNU_EPJ_SLIPIDS/em.tpr',\n",
" './SIMULATIONS/NACHRA7_PNU_EPJ_SLIPIDS/em.gro')\n",
"check_universe_ring_penetration(universe)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"In frame 0\n",
"%lipid ring pentration\n"
]
}
],
"source": [
"universe = mda.Universe('./SIMULATIONS/NACHRA7_PNU_EPJ_SLIPIDS/alchem.tpr',\n",
" './SIMULATIONS/NACHRA7_PNU_EPJ_SLIPIDS/alchem.gro')\n",
"check_universe_ring_penetration(universe)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "mda_py38",
"language": "python",
"name": "jupyter"
},
"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.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment