Skip to content

Instantly share code, notes, and snippets.

@zhang-ivy
Last active August 6, 2020 17:12
Show Gist options
  • Save zhang-ivy/1dba04332bd619e7d29d018804ae2274 to your computer and use it in GitHub Desktop.
Save zhang-ivy/1dba04332bd619e7d29d018804ae2274 to your computer and use it in GitHub Desktop.
prep_ncmc_movie
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get small molecule sdf"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from openeye import oechem\n",
"\n",
"ifs = oechem.oemolistream(\"inhibitor.pdb\")\n",
"ofs = oechem.oemolostream(\"inhibitor.sdf\")\n",
"\n",
"ifs.SetFormat(oechem.OEFormat_PDB)\n",
"ofs.SetFormat(oechem.OEFormat_SDF)\n",
"\n",
"for mol in ifs.GetOEGraphMols():\n",
" oechem.OEWriteMolecule(ofs, mol)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate htf"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Set up logger\n",
"import logging\n",
"_logger = logging.getLogger()\n",
"_logger.setLevel(logging.DEBUG)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:utils.openeye:molecule does not have unique atom names. Generating now...\n",
"INFO:utils.openeye:molecule has unique atom names already\n",
"INFO:root:solvating at 0.15 M using tip3p\n",
"INFO:root:solvating at 0.15 M using tip3p\n",
"INFO:openmmforcefields.generators.template_generators:Requested to generate parameters for residue <Residue 297 () of chain 2>\n",
"INFO:openmmforcefields.generators.template_generators:Generating a residue template for [H]c1c(c(c(c(c1F)[H])[C@]2(C(C(C([N@@]2C3=NC4=C(C(=NN4C(=C3[H])[H])[H])N([H])C(=O)N5C([C@@](C(C5([H])[H])([H])[H])([H])O[H])([H])[H])([H])[H])([H])[H])([H])[H])[H])F)[H]\n",
"INFO:proposal_generator:\tConducting polymer point mutation proposal...\n",
"INFO:proposal_generator:Using matching_criterion to chose best atom map\n",
"INFO:proposal_generator:Scaffold has symmetry of 0\n",
"WARNING:proposal_generator:Two molecules are not similar to have a common scaffold\n",
"WARNING:proposal_generator:Proceeding with direct mapping of molecules, but please check atom mapping and the geometry of the ligands.\n",
"INFO:proposal_generator:len [{}, {}, {8: 6}, {9: 6}, {10: 6}]\n",
"INFO:proposal_generator:{}\n",
"INFO:proposal_generator:{}\n",
"INFO:proposal_generator:{8: 6}\n",
"INFO:proposal_generator:{9: 6}\n",
"INFO:proposal_generator:{10: 6}\n",
"INFO:proposal_generator:Returning map that best satisfies matching_criterion\n",
"INFO:proposal_generator:Finding best map using matching_criterion name\n",
"INFO:proposal_generator:{}\n",
"INFO:geometry:propose: performing forward proposal\n",
"INFO:geometry:propose: unique new atoms detected; proceeding to _logp_propose...\n",
"INFO:geometry:Conducting forward proposal...\n",
"INFO:geometry:Computing proposal order with NetworkX...\n",
"INFO:geometry:number of atoms to be placed: 5\n",
"INFO:geometry:Atom index proposal order is [2611, 2614, 2613, 2612, 2615]\n",
"INFO:geometry:omitted_bonds: []\n",
"INFO:geometry:direction of proposal is forward; creating atoms_with_positions and new positions from old system/topology...\n",
"INFO:geometry:creating growth system...\n",
"INFO:geometry:\tcreating bond force...\n",
"INFO:geometry:\tthere are 2423 bonds in reference force.\n",
"INFO:geometry:\tcreating angle force...\n",
"INFO:geometry:\tthere are 8641 angles in reference force.\n",
"INFO:geometry:\tcreating torsion force...\n",
"INFO:geometry:\tcreating extra torsions force...\n",
"INFO:geometry:\tthere are 16054 torsions in reference force.\n",
"INFO:geometry:\tcreating nonbonded force...\n",
"INFO:geometry:\t\tgrabbing reference nonbonded method, cutoff, switching function, switching distance...\n",
"INFO:geometry:\t\tcreating nonbonded exception force (i.e. custom bond for 1,4s)...\n",
"INFO:geometry:\t\tlooping through exceptions calculating growth indices, and adding appropriate interactions to custom bond force.\n",
"INFO:geometry:\t\tthere are 73489 in the reference Nonbonded force\n",
"WARNING:geometry:\t\t\t\t\tchiral atom <Atom 2606 (CA) of chain 0 residue 166 (CYS)> with neighbors [<Atom 2607 (C) of chain 0 residue 166 (CYS)>, <Atom 2605 (N) of chain 0 residue 166 (CYS)>, <Atom 2610 (HA) of chain 0 residue 166 (CYS)>, <Atom 2611 (CB) of chain 0 residue 166 (CYS)>] is surrounded by 3 core neighbors. omitting chirality bias torsion\n",
"INFO:geometry:Neglected angle terms : []\n",
"INFO:geometry:omitted_growth_terms: {'bonds': [], 'angles': [], 'torsions': [], '1,4s': []}\n",
"INFO:geometry:extra torsions: {}\n",
"INFO:geometry:neglected angle terms include []\n",
"INFO:geometry:log probability choice of torsions and atom order: -6.4457198193855785\n",
"INFO:geometry:creating platform, integrators, and contexts; setting growth parameter\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:setting atoms_with_positions context new positions\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:There are 5 new atoms\n",
"INFO:geometry:\treduced angle potential = 0.6793831317322927.\n",
"INFO:geometry:\treduced angle potential = 0.7074218930654235.\n",
"INFO:geometry:\treduced angle potential = 0.03453944289459454.\n",
"INFO:geometry:\treduced angle potential = 0.06148813429814449.\n",
"INFO:geometry:\treduced angle potential = 1.4016947314305477.\n",
"INFO:geometry:\tbeginning construction of no_nonbonded final system...\n",
"INFO:geometry:\tinitial no-nonbonded final system forces ['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat']\n",
"INFO:geometry:\tfinal no-nonbonded final system forces dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce'])\n",
"INFO:geometry:\tthere are 2423 bond forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 8641 angle forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 16054 torsion forces in the no-nonbonded final system\n",
"INFO:geometry:forward final system defined with 0 neglected angles.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced potential before atom placement: 12119.007527980379\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced energy added from growth system: 17.1324098013441\n",
"INFO:geometry:final reduced energy 12136.139936713791\n",
"INFO:geometry:sum of energies: 12136.139937781723\n",
"INFO:geometry:magnitude of difference in the energies: 1.067931723497395e-06\n",
"INFO:geometry:Final logp_proposal: 31.257119159893886\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"added energy components: [('CustomBondForce', 0.1788297299321343), ('CustomAngleForce', 7.0777452275938275), ('CustomTorsionForce', 9.317535099123424), ('CustomBondForce', 0.5582997446947136)]\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:logp_reverse: performing reverse proposal\n",
"INFO:geometry:logp_reverse: unique new atoms detected; proceeding to _logp_propose...\n",
"INFO:geometry:Conducting forward proposal...\n",
"INFO:geometry:Computing proposal order with NetworkX...\n",
"INFO:geometry:number of atoms to be placed: 1\n",
"INFO:geometry:Atom index proposal order is [2610]\n",
"INFO:geometry:omitted_bonds: []\n",
"INFO:geometry:direction of proposal is reverse; creating atoms_with_positions from old system/topology\n",
"INFO:geometry:creating growth system...\n",
"INFO:geometry:\tcreating bond force...\n",
"INFO:geometry:\tthere are 2421 bonds in reference force.\n",
"INFO:geometry:\tcreating angle force...\n",
"INFO:geometry:\tthere are 8634 angles in reference force.\n",
"INFO:geometry:\tcreating torsion force...\n",
"INFO:geometry:\tcreating extra torsions force...\n",
"INFO:geometry:\tthere are 16029 torsions in reference force.\n",
"INFO:geometry:\tcreating nonbonded force...\n",
"INFO:geometry:\t\tgrabbing reference nonbonded method, cutoff, switching function, switching distance...\n",
"INFO:geometry:\t\tcreating nonbonded exception force (i.e. custom bond for 1,4s)...\n",
"INFO:geometry:\t\tlooping through exceptions calculating growth indices, and adding appropriate interactions to custom bond force.\n",
"INFO:geometry:\t\tthere are 73466 in the reference Nonbonded force\n",
"INFO:geometry:Neglected angle terms : []\n",
"INFO:geometry:omitted_growth_terms: {'bonds': [], 'angles': [], 'torsions': [], '1,4s': []}\n",
"INFO:geometry:extra torsions: {}\n",
"INFO:geometry:neglected angle terms include []\n",
"INFO:geometry:log probability choice of torsions and atom order: -0.6931471805599453\n",
"INFO:geometry:creating platform, integrators, and contexts; setting growth parameter\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:setting atoms_with_positions context old positions\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:There are 1 new atoms\n",
"INFO:geometry:\treduced angle potential = 1.0299774016188954.\n",
"INFO:geometry:\tbeginning construction of no_nonbonded final system...\n",
"INFO:geometry:\tinitial no-nonbonded final system forces ['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat']\n",
"INFO:geometry:\tfinal no-nonbonded final system forces dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce'])\n",
"INFO:geometry:\tthere are 2421 bond forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 8634 angle forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 16029 torsion forces in the no-nonbonded final system\n",
"INFO:geometry:reverse final system defined with 0 neglected angles.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced potential before atom placement: 12119.60516216377\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced energy added from growth system: -0.07650105018990651\n",
"INFO:geometry:final reduced energy 12119.52866015684\n",
"INFO:geometry:sum of energies: 12119.528661113582\n",
"INFO:geometry:magnitude of difference in the energies: 9.567403343874181e-07\n",
"INFO:geometry:Final logp_proposal: 6.068315745841434\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"added energy components: [('CustomBondForce', 0.0), ('CustomAngleForce', 1.1923846415675037), ('CustomTorsionForce', 1.1151859981354997), ('CustomBondForce', -2.384071689892909)]\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:relative:Beginning nonbonded method, total particle, barostat, and exceptions retrieval...\n",
"INFO:relative:Old system forces: dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat'])\n",
"INFO:relative:New system forces: dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat'])\n",
"INFO:relative:No unknown forces.\n",
"INFO:relative:Nonbonded method to be used (i.e. from old system): 4\n",
"INFO:relative:Adding and mapping old atoms to hybrid system...\n",
"INFO:relative:Adding and mapping new atoms to hybrid system...\n",
"INFO:relative:Added MonteCarloBarostat.\n",
"INFO:relative:getDefaultPeriodicBoxVectors added to hybrid: [Quantity(value=Vec3(x=8.221, y=0.0, z=0.0), unit=nanometer), Quantity(value=Vec3(x=0.0, y=8.221, z=0.0), unit=nanometer), Quantity(value=Vec3(x=0.0, y=0.0, z=8.221), unit=nanometer)]\n",
"INFO:relative:Determined atom classes.\n",
"INFO:relative:Generating old system exceptions dict...\n",
"INFO:relative:Generating new system exceptions dict...\n",
"INFO:relative:Handling constraints...\n",
"INFO:relative:Handling virtual sites...\n",
"INFO:relative:\t_handle_virtual_sites: numVirtualSites: 0\n",
"INFO:relative:Adding bond force terms...\n",
"INFO:relative:Adding angle force terms...\n",
"INFO:relative:Adding torsion force terms...\n",
"INFO:relative:Adding nonbonded force terms...\n",
"INFO:relative:\t_add_nonbonded_force_terms: <simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x2ad9a10e4600> > added to hybrid system\n",
"INFO:relative:\t_add_nonbonded_force_terms: nonbonded_method is PME or Ewald\n",
"INFO:relative:\t_add_nonbonded_force_terms: 4 added to standard nonbonded force\n",
"INFO:relative:\t_add_nonbonded_force_terms: 2 added to sterics_custom_nonbonded force\n",
"INFO:relative:\t_add_nonbonded_force_terms: <simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x2ad9a10e43f0> > added to hybrid system\n",
"INFO:relative:Handling harmonic bonds...\n",
"INFO:relative:\thandle_harmonic_bonds: looping through old_system to add relevant terms...\n",
"INFO:relative:\thandle_harmonic_bonds: looping through new_system to add relevant terms...\n",
"INFO:relative:Handling harmonic angles...\n",
"INFO:relative:\thandle_harmonic_angles: looping through old_system to add relevant terms...\n",
"INFO:relative:\thandle_harmonic_angles: looping through new_system to add relevant terms...\n",
"INFO:relative:Handling torsion forces...\n",
"INFO:relative:\thandle_periodic_torsion_forces: looping through old_system to add relevant terms...\n",
"INFO:relative:\thandle_periodic_torsion_forces: looping through new_system to add relevant terms...\n",
"INFO:relative:Handling nonbonded forces...\n",
"INFO:relative:\thandle_nonbonded: looping through all particles in hybrid...\n",
"INFO:relative:\thandle_nonbonded: Handling Interaction Groups...\n",
"INFO:relative:\thandle_nonbonded: Handling Hybrid Exceptions...\n",
"INFO:relative:\thandle_nonbonded: Handling Original Exceptions...\n",
"INFO:relative:Handling unique_new/old interaction exceptions...\n",
"INFO:relative:There are old or new system exceptions...proceeding.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tHarmonicBondForce: 1586.6864162448433\n",
"\t\t\tHarmonicAngleForce: 4286.78466551953\n",
"\t\t\tPeriodicTorsionForce: 6248.441650659238\n",
"\t\t\tNonbondedForce: -176969.5206284497\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tCustomBondForce: 8.597306569829044\n",
"\t\t\tHarmonicBondForce: 1578.267939404946\n",
"\t\t\tCustomAngleForce: 18.922580479241798\n",
"\t\t\tHarmonicAngleForce: 4274.939830267883\n",
"\t\t\tCustomTorsionForce: 9.425421602998034\n",
"\t\t\tPeriodicTorsionForce: 6248.333785959808\n",
"\t\t\tNonbondedForce: -176907.95420332515\n",
"\t\t\tCustomNonbondedForce: -61.00971938570378\n",
"\t\t\tCustomBondForce: 0.0\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tCustomBondForce: 8.597306569829044\n",
"\t\t\tHarmonicBondForce: 1578.267939404946\n",
"\t\t\tCustomAngleForce: 18.922580479241798\n",
"\t\t\tHarmonicAngleForce: 4274.939830267883\n",
"\t\t\tCustomTorsionForce: 8.827786766234517\n",
"\t\t\tPeriodicTorsionForce: 6248.333785959808\n",
"\t\t\tNonbondedForce: -176899.77889609477\n",
"\t\t\tCustomNonbondedForce: 2015.3988151315605\n",
"\t\t\tCustomBondForce: 0.0\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tHarmonicBondForce: 1586.8652459747752\n",
"\t\t\tHarmonicAngleForce: 4292.670026105557\n",
"\t\t\tPeriodicTorsionForce: 6256.046365131089\n",
"\t\t\tNonbondedForce: -174881.9944453129\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:proposal_generator:\tConducting polymer point mutation proposal...\n",
"INFO:proposal_generator:Using matching_criterion to chose best atom map\n",
"INFO:proposal_generator:Scaffold has symmetry of 0\n",
"WARNING:proposal_generator:Two molecules are not similar to have a common scaffold\n",
"WARNING:proposal_generator:Proceeding with direct mapping of molecules, but please check atom mapping and the geometry of the ligands.\n",
"INFO:proposal_generator:len [{}, {}, {8: 6}, {9: 6}, {10: 6}]\n",
"INFO:proposal_generator:{}\n",
"INFO:proposal_generator:{}\n",
"INFO:proposal_generator:{8: 6}\n",
"INFO:proposal_generator:{9: 6}\n",
"INFO:proposal_generator:{10: 6}\n",
"INFO:proposal_generator:Returning map that best satisfies matching_criterion\n",
"INFO:proposal_generator:Finding best map using matching_criterion name\n",
"INFO:proposal_generator:{}\n",
"INFO:geometry:propose: performing forward proposal\n",
"INFO:geometry:propose: unique new atoms detected; proceeding to _logp_propose...\n",
"INFO:geometry:Conducting forward proposal...\n",
"INFO:geometry:Computing proposal order with NetworkX...\n",
"INFO:geometry:number of atoms to be placed: 5\n",
"INFO:geometry:Atom index proposal order is [2611, 2614, 2612, 2615, 2613]\n",
"INFO:geometry:omitted_bonds: []\n",
"INFO:geometry:direction of proposal is forward; creating atoms_with_positions and new positions from old system/topology...\n",
"INFO:geometry:creating growth system...\n",
"INFO:geometry:\tcreating bond force...\n",
"INFO:geometry:\tthere are 2458 bonds in reference force.\n",
"INFO:geometry:\tcreating angle force...\n",
"INFO:geometry:\tthere are 8743 angles in reference force.\n",
"INFO:geometry:\tcreating torsion force...\n",
"INFO:geometry:\tcreating extra torsions force...\n",
"INFO:geometry:\tthere are 16213 torsions in reference force.\n",
"INFO:geometry:\tcreating nonbonded force...\n",
"INFO:geometry:\t\tgrabbing reference nonbonded method, cutoff, switching function, switching distance...\n",
"INFO:geometry:\t\tcreating nonbonded exception force (i.e. custom bond for 1,4s)...\n",
"INFO:geometry:\t\tlooping through exceptions calculating growth indices, and adding appropriate interactions to custom bond force.\n",
"INFO:geometry:\t\tthere are 73861 in the reference Nonbonded force\n",
"WARNING:geometry:\t\t\t\t\tchiral atom <Atom 2606 (CA) of chain 0 residue 166 (CYS)> with neighbors [<Atom 2607 (C) of chain 0 residue 166 (CYS)>, <Atom 2605 (N) of chain 0 residue 166 (CYS)>, <Atom 2610 (HA) of chain 0 residue 166 (CYS)>, <Atom 2611 (CB) of chain 0 residue 166 (CYS)>] is surrounded by 3 core neighbors. omitting chirality bias torsion\n",
"INFO:geometry:Neglected angle terms : []\n",
"INFO:geometry:omitted_growth_terms: {'bonds': [], 'angles': [], 'torsions': [], '1,4s': []}\n",
"INFO:geometry:extra torsions: {}\n",
"INFO:geometry:neglected angle terms include []\n",
"INFO:geometry:log probability choice of torsions and atom order: -6.733401891837359\n",
"INFO:geometry:creating platform, integrators, and contexts; setting growth parameter\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:setting atoms_with_positions context new positions\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:There are 5 new atoms\n",
"INFO:geometry:\treduced angle potential = 1.0162203615594578.\n",
"INFO:geometry:\treduced angle potential = 1.0554286619518236.\n",
"INFO:geometry:\treduced angle potential = 1.6225427964638328.\n",
"INFO:geometry:\treduced angle potential = 0.07391339712277141.\n",
"INFO:geometry:\treduced angle potential = 0.07220977573850834.\n",
"INFO:geometry:\tbeginning construction of no_nonbonded final system...\n",
"INFO:geometry:\tinitial no-nonbonded final system forces ['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat']\n",
"INFO:geometry:\tfinal no-nonbonded final system forces dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce'])\n",
"INFO:geometry:\tthere are 2458 bond forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 8743 angle forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 16213 torsion forces in the no-nonbonded final system\n",
"INFO:geometry:forward final system defined with 0 neglected angles.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced potential before atom placement: 12283.325417704327\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced energy added from growth system: 24.8449619707893\n",
"INFO:geometry:final reduced energy 12308.170377137325\n",
"INFO:geometry:sum of energies: 12308.170379675117\n",
"INFO:geometry:magnitude of difference in the energies: 2.537791587542415e-06\n",
"INFO:geometry:Final logp_proposal: 31.2529965071565\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"added energy components: [('CustomBondForce', 0.08001307083457838), ('CustomAngleForce', 5.259671402486339), ('CustomTorsionForce', 11.23827319051712), ('CustomBondForce', 8.267004306951257)]\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:logp_reverse: performing reverse proposal\n",
"INFO:geometry:logp_reverse: unique new atoms detected; proceeding to _logp_propose...\n",
"INFO:geometry:Conducting forward proposal...\n",
"INFO:geometry:Computing proposal order with NetworkX...\n",
"INFO:geometry:number of atoms to be placed: 1\n",
"INFO:geometry:Atom index proposal order is [2610]\n",
"INFO:geometry:omitted_bonds: []\n",
"INFO:geometry:direction of proposal is reverse; creating atoms_with_positions from old system/topology\n",
"INFO:geometry:creating growth system...\n",
"INFO:geometry:\tcreating bond force...\n",
"INFO:geometry:\tthere are 2456 bonds in reference force.\n",
"INFO:geometry:\tcreating angle force...\n",
"INFO:geometry:\tthere are 8736 angles in reference force.\n",
"INFO:geometry:\tcreating torsion force...\n",
"INFO:geometry:\tcreating extra torsions force...\n",
"INFO:geometry:\tthere are 16188 torsions in reference force.\n",
"INFO:geometry:\tcreating nonbonded force...\n",
"INFO:geometry:\t\tgrabbing reference nonbonded method, cutoff, switching function, switching distance...\n",
"INFO:geometry:\t\tcreating nonbonded exception force (i.e. custom bond for 1,4s)...\n",
"INFO:geometry:\t\tlooping through exceptions calculating growth indices, and adding appropriate interactions to custom bond force.\n",
"INFO:geometry:\t\tthere are 73838 in the reference Nonbonded force\n",
"INFO:geometry:Neglected angle terms : []\n",
"INFO:geometry:omitted_growth_terms: {'bonds': [], 'angles': [], 'torsions': [], '1,4s': []}\n",
"INFO:geometry:extra torsions: {}\n",
"INFO:geometry:neglected angle terms include []\n",
"INFO:geometry:log probability choice of torsions and atom order: -0.6931471805599453\n",
"INFO:geometry:creating platform, integrators, and contexts; setting growth parameter\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:setting atoms_with_positions context old positions\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:There are 1 new atoms\n",
"INFO:geometry:\treduced angle potential = 0.017863316105502774.\n",
"INFO:geometry:\tbeginning construction of no_nonbonded final system...\n",
"INFO:geometry:\tinitial no-nonbonded final system forces ['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat']\n",
"INFO:geometry:\tfinal no-nonbonded final system forces dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce'])\n",
"INFO:geometry:\tthere are 2456 bond forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 8736 angle forces in the no-nonbonded final system\n",
"INFO:geometry:\tthere are 16188 torsion forces in the no-nonbonded final system\n",
"INFO:geometry:reverse final system defined with 0 neglected angles.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced potential before atom placement: 12283.921842317495\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:geometry:total reduced energy added from growth system: -0.07650105018990651\n",
"INFO:geometry:final reduced energy 12283.845340310561\n",
"INFO:geometry:sum of energies: 12283.845341267306\n",
"INFO:geometry:magnitude of difference in the energies: 9.567439723662252e-07\n",
"INFO:geometry:Final logp_proposal: 5.618041238207212\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"added energy components: [('CustomBondForce', 0.0), ('CustomAngleForce', 1.1923846415675037), ('CustomTorsionForce', 1.1151859981354997), ('CustomBondForce', -2.384071689892909)]\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO:relative:Beginning nonbonded method, total particle, barostat, and exceptions retrieval...\n",
"INFO:relative:Old system forces: dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat'])\n",
"INFO:relative:New system forces: dict_keys(['HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat'])\n",
"INFO:relative:No unknown forces.\n",
"INFO:relative:Nonbonded method to be used (i.e. from old system): 4\n",
"INFO:relative:Adding and mapping old atoms to hybrid system...\n",
"INFO:relative:Adding and mapping new atoms to hybrid system...\n",
"INFO:relative:Added MonteCarloBarostat.\n",
"INFO:relative:getDefaultPeriodicBoxVectors added to hybrid: [Quantity(value=Vec3(x=8.221, y=0.0, z=0.0), unit=nanometer), Quantity(value=Vec3(x=0.0, y=8.221, z=0.0), unit=nanometer), Quantity(value=Vec3(x=0.0, y=0.0, z=8.221), unit=nanometer)]\n",
"INFO:relative:Determined atom classes.\n",
"INFO:relative:Generating old system exceptions dict...\n",
"INFO:relative:Generating new system exceptions dict...\n",
"INFO:relative:Handling constraints...\n",
"INFO:relative:Handling virtual sites...\n",
"INFO:relative:\t_handle_virtual_sites: numVirtualSites: 0\n",
"INFO:relative:Adding bond force terms...\n",
"INFO:relative:Adding angle force terms...\n",
"INFO:relative:Adding torsion force terms...\n",
"INFO:relative:Adding nonbonded force terms...\n",
"INFO:relative:\t_add_nonbonded_force_terms: <simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x2ad9e3a43420> > added to hybrid system\n",
"INFO:relative:\t_add_nonbonded_force_terms: nonbonded_method is PME or Ewald\n",
"INFO:relative:\t_add_nonbonded_force_terms: 4 added to standard nonbonded force\n",
"INFO:relative:\t_add_nonbonded_force_terms: 2 added to sterics_custom_nonbonded force\n",
"INFO:relative:\t_add_nonbonded_force_terms: <simtk.openmm.openmm.CustomNonbondedForce; proxy of <Swig Object of type 'OpenMM::CustomNonbondedForce *' at 0x2ad9e3a43090> > added to hybrid system\n",
"INFO:relative:Handling harmonic bonds...\n",
"INFO:relative:\thandle_harmonic_bonds: looping through old_system to add relevant terms...\n",
"INFO:relative:\thandle_harmonic_bonds: looping through new_system to add relevant terms...\n",
"INFO:relative:Handling harmonic angles...\n",
"INFO:relative:\thandle_harmonic_angles: looping through old_system to add relevant terms...\n",
"INFO:relative:\thandle_harmonic_angles: looping through new_system to add relevant terms...\n",
"INFO:relative:Handling torsion forces...\n",
"INFO:relative:\thandle_periodic_torsion_forces: looping through old_system to add relevant terms...\n",
"INFO:relative:\thandle_periodic_torsion_forces: looping through new_system to add relevant terms...\n",
"INFO:relative:Handling nonbonded forces...\n",
"INFO:relative:\thandle_nonbonded: looping through all particles in hybrid...\n",
"INFO:relative:\thandle_nonbonded: Handling Interaction Groups...\n",
"INFO:relative:\thandle_nonbonded: Handling Hybrid Exceptions...\n",
"INFO:relative:\thandle_nonbonded: Handling Original Exceptions...\n",
"INFO:relative:Handling unique_new/old interaction exceptions...\n",
"INFO:relative:There are old or new system exceptions...proceeding.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"conducting subsequent work with the following platform: CUDA\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tHarmonicBondForce: 1608.1379485146708\n",
"\t\t\tHarmonicAngleForce: 4397.414686393739\n",
"\t\t\tPeriodicTorsionForce: 6280.676777668922\n",
"\t\t\tNonbondedForce: -186046.63272028015\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tCustomBondForce: 8.597306569829044\n",
"\t\t\tHarmonicBondForce: 1599.6206550156762\n",
"\t\t\tCustomAngleForce: 18.922580479241798\n",
"\t\t\tHarmonicAngleForce: 4383.751777316984\n",
"\t\t\tCustomTorsionForce: 9.491092511012464\n",
"\t\t\tPeriodicTorsionForce: 6282.423978799392\n",
"\t\t\tNonbondedForce: -185969.93943816153\n",
"\t\t\tCustomNonbondedForce: -68.42786941568974\n",
"\t\t\tCustomBondForce: 0.0\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tCustomBondForce: 8.597306569829044\n",
"\t\t\tHarmonicBondForce: 1599.6206550156762\n",
"\t\t\tCustomAngleForce: 18.922580479241798\n",
"\t\t\tHarmonicAngleForce: 4383.751777316984\n",
"\t\t\tCustomTorsionForce: 8.894667329726287\n",
"\t\t\tPeriodicTorsionForce: 6282.423978799392\n",
"\t\t\tNonbondedForce: -185957.77056570715\n",
"\t\t\tCustomNonbondedForce: 8513.98059371616\n",
"\t\t\tCustomBondForce: 0.0\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n",
"conducting subsequent work with the following platform: CUDA\n",
"\t\t\tHarmonicBondForce: 1608.2179615855052\n",
"\t\t\tHarmonicAngleForce: 4401.481973154659\n",
"\t\t\tPeriodicTorsionForce: 6290.20343841989\n",
"\t\t\tNonbondedForce: -177441.40495286722\n",
"\t\t\tAndersenThermostat: 0.0\n",
"\t\t\tMonteCarloBarostat: 0.0\n"
]
}
],
"source": [
"# Solvent\n",
"from perses.app.relative_point_mutation_setup import PointMutationExecutor\n",
"from simtk import unit\n",
"solvent_delivery = PointMutationExecutor(\"kinase.pdb\", \n",
" '1', # First chain is the barstar one\n",
" '667', \n",
" 'CYS',\n",
" ligand_file=\"inhibitor.sdf\",\n",
" ionic_strength=0.15*unit.molar\n",
" )\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Check maps"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def render_protein_residue_atom_mapping(topology_proposal, filename):\n",
" \"\"\"\n",
" wrap the `render_atom_mapping` method around protein point mutation topologies.\n",
" TODO : make modification to `render_atom_mapping` so that the backbone atoms are not written in the output.\n",
" \n",
" arguments\n",
" topology_proposal : perses.rjmc.topology_proposal.TopologyProposal object\n",
" topology proposal of protein mutation\n",
" filename : str\n",
" filename to write the map\n",
" \"\"\"\n",
" from perses.utils.smallmolecules import render_atom_mapping\n",
" oe_res_maps = {}\n",
" for omm_new_idx, omm_old_idx in topology_proposal._new_to_old_atom_map.items():\n",
" if omm_new_idx in topology_proposal._new_topology.residue_to_oemol_map.keys():\n",
" try:\n",
" oe_res_maps[topology_proposal._new_topology.residue_to_oemol_map[omm_new_idx]] = topology_proposal._old_topology.residue_to_oemol_map[omm_old_idx]\n",
" except:\n",
" pass\n",
" \n",
" render_atom_mapping(filename, topology_proposal._old_topology.residue_oemol, topology_proposal._new_topology.residue_oemol, oe_res_maps)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"render_protein_residue_atom_mapping(solvent_delivery.get_apo_htf()._topology_proposal, \"apo_map.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Save htfs as pickles"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"import os\n",
"out_dir = '/data/chodera/zhangi/ncmc_movie/0/'\n",
"pickle.dump(solvent_delivery.get_apo_htf(), open(os.path.join(out_dir, \"0_apo.pickle\"), \"wb\" ))\n",
"pickle.dump(solvent_delivery.get_complex_htf(), open(os.path.join(out_dir, \"0_complex.pickle\"), \"wb\" ))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Save old and new positions for apo and complex to check geometry"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"import mdtraj as md\n",
"import numpy as np\n",
"htfs_t = [solvent_delivery.get_apo_htf(), solvent_delivery.get_complex_htf()]\n",
"\n",
"top_old = md.Topology.from_openmm(htfs_t[0]._topology_proposal.old_topology)\n",
"top_new = md.Topology.from_openmm(htfs_t[0]._topology_proposal.new_topology)\n",
"traj = md.Trajectory(np.array(htfs_t[0].old_positions(htfs_t[0].hybrid_positions)), top_old)\n",
"traj.save(\"0_apo_old.pdb\")\n",
"traj = md.Trajectory(np.array(htfs_t[0].new_positions(htfs_t[0].hybrid_positions)), top_new)\n",
"traj.save(\"0_apo_new.pdb\")\n",
"\n",
"top_old = md.Topology.from_openmm(htfs_t[1]._topology_proposal.old_topology)\n",
"top_new = md.Topology.from_openmm(htfs_t[1]._topology_proposal.new_topology)\n",
"traj = md.Trajectory(np.array(htfs_t[1].old_positions(htfs_t[1].hybrid_positions)), top_old)\n",
"traj.save(\"0_complex_old.pdb\")\n",
"traj = md.Trajectory(np.array(htfs_t[1].new_positions(htfs_t[1].hybrid_positions)), top_new)\n",
"traj.save(\"0_complex_new.pdb\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment