updated mapping code and logger after recent bugfixes ***NO_CI***

Author: jmichel80

nodes/playground/prepareFEP.ipynb

@@ -38,10 +38,22 @@
     "    \"\"\" Human readable report on atoms used for the mapping.\"\"\"\n",
     "    atoms_in_A = list(mapping.keys())\n",
     "    stream = open('somd.mapping','w')\n",
+    "    atAdone = []\n",
+    "    atBdone= []\n",
     "    for atAidx in atoms_in_A:\n",
     "        atA = ligA._sire_molecule.select(atAidx)\n",
     "        atB = ligB._sire_molecule.select(mapping[atAidx])\n",
-    "        stream.write(\"%s --> %s\\n\" % (atA.name(),atB.name()))\n",
+    "        stream.write(\"%s %s --> %s %s\\n\" % (atA.index(), atA.name(),atB.index(), atB.name()))\n",
+    "        atAdone.append(atA)\n",
+    "        atBdone.append(atB)\n",
+    "    for atom in ligA._sire_molecule.atoms():\n",
+    "        if atom in atAdone:\n",
+    "            continue\n",
+    "        stream.write(\"%s %s --> dummy\\n\" % (atom.index(), atom.name()))\n",
+    "    for atom in ligB._sire_molecule.atoms():\n",
+    "        if atom in atBdone:\n",
+    "            continue\n",
+    "        stream.write(\"dummy --> %s %s\\n\" % (atom.index(), atom.name()))\n",
     "    stream.close()"
    ]
   },
@@ -72,7 +84,7 @@
    "source": [
     "node.addInput(\"input1\", BSS.Gateway.FileSet(help=\"A topology and coordinates file\"))\n",
     "node.addInput(\"input2\", BSS.Gateway.FileSet(help=\"A topology and coordinates file\"))\n",
-    "node.addInput(\"prematch\", BSS.Gateway.String(help=\"list of atom indices that are matched between input1 and input2. Syntax is of the format 1-3,4-8,9-11...\", default=\"\"))\n",
+    "node.addInput(\"prematch\", BSS.Gateway.String(help=\"list of atom indices that are matched between input2 and input1. Syntax is of the format 1-3,4-8,9-11...\", default=\"\"))\n",
     "node.addInput(\"output\", BSS.Gateway.String(help=\"The root name for the files describing the perturbation input1->input2.\"))"
    ]
   },
@@ -103,11 +115,12 @@
     "# Optional input, dictionary of Atom indices that should be matched in the search. \n",
     "prematch = {}\n",
     "prematchstring = node.getInput(\"prematch\")\n",
-    "entries = prematchstring.split(\",\")\n",
-    "for entry in entries:\n",
-    "    idxA, idxB = entry.split(\"-\")\n",
-    "    prematch[ AtomIdx( int(idxA)) ] = AtomIdx( int(idxB) )\n",
-    "print (prematch)"
+    "if len(prematchstring) > 0: \n",
+    "    entries = prematchstring.split(\",\")\n",
+    "    for entry in entries:\n",
+    "        idxA, idxB = entry.split(\"-\")\n",
+    "        prematch[ AtomIdx( int(idxA)) ] = AtomIdx( int(idxB) )\n",
+    "#print (prematch)"
    ]
   },
   {
@@ -148,22 +161,11 @@
    "outputs": [],
    "source": [
     "# Return a maximum of 10 matches, scored by RMSD and sorted from best to worst.\n",
-    "mappings = BSS.Align.matchAtoms(lig1, lig2, matches=10, prematch=prematch)\n",
+    "mappings = BSS.Align.matchAtoms(lig1, lig2, matches=10, prematch=prematch, scoring_function=\"RMSD\", timeout=1*BSS.Units.Time.minute)\n",
     "# We retain the top mapping\n",
     "mapping = mappings[0]\n",
-    "print (len(mappings))\n",
-    "print (mappings)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "for x in range(0,len(mappings)):\n",
-    "    print (mappings[x])\n",
-    "    print (\"-----------\")"
+    "#print (len(mappings))\n",
+    "#print (mappings)"
    ]
   },
   {
@@ -267,7 +269,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.5.5"
+   "version": "3.7.1"
   }
  },
  "nbformat": 4,

nodes/playground/prepareFEP.py

@@ -1,4 +1,4 @@
-
+#!/usr/bin/env python
 # coding: utf-8
 
 # Author: Julien Michel
@@ -8,72 +8,86 @@
 # # PrepareFEP
 # Loads a pair of input files, perform mapping between the first molecule of each input. Write down input files for a SOMD FEP calculation.
 
-# In[8]:
+# In[1]:
 
 
 import BioSimSpace as BSS
 import os
 from Sire.Mol import AtomIdx
 
 
-# In[9]:
+# In[2]:
 
 
 def writeLog(ligA, ligB, mapping):
     """ Human readable report on atoms used for the mapping."""
     atoms_in_A = list(mapping.keys())
     stream = open('somd.mapping','w')
+    atAdone = []
+    atBdone= []
     for atAidx in atoms_in_A:
         atA = ligA._sire_molecule.select(atAidx)
         atB = ligB._sire_molecule.select(mapping[atAidx])
-        stream.write("%s --> %s\n" % (atA.name(),atB.name()))
+        stream.write("%s %s --> %s %s\n" % (atA.index(), atA.name(),atB.index(), atB.name()))
+        atAdone.append(atA)
+        atBdone.append(atB)
+    for atom in ligA._sire_molecule.atoms():
+        if atom in atAdone:
+            continue
+        stream.write("%s %s --> dummy\n" % (atom.index(), atom.name()))
+    for atom in ligB._sire_molecule.atoms():
+        if atom in atBdone:
+            continue
+        stream.write("dummy --> %s %s\n" % (atom.index(), atom.name()))
     stream.close()
 
 
-# In[10]:
+# In[3]:
 
 
 node = BSS.Gateway.Node("A node to generate input files for a SOMD relative free energy calculation.")
 
 
-# In[11]:
+# In[4]:
 
 
 node.addAuthor(name="Julien Michel", email="julien.michel@ed.ac.uk", affiliation="University of Edinburgh")
 node.setLicense("GPLv3")
 
 
-# In[12]:
+# In[5]:
 
 
 node.addInput("input1", BSS.Gateway.FileSet(help="A topology and coordinates file"))
 node.addInput("input2", BSS.Gateway.FileSet(help="A topology and coordinates file"))
-node.addInput("prematch", BSS.Gateway.String(help="list of atom indices that are matched between input1 and input2. Syntax is of the format 1-3,4-8,9-11...", default=""))
+node.addInput("prematch", BSS.Gateway.String(help="list of atom indices that are matched between input2 and input1. Syntax is of the format 1-3,4-8,9-11...", default=""))
 node.addInput("output", BSS.Gateway.String(help="The root name for the files describing the perturbation input1->input2."))
 
 
-# In[13]:
+# In[6]:
 
 
 node.addOutput("nodeoutput", BSS.Gateway.FileSet(help="SOMD input files for a perturbation of input1->input2."))
 
 
-# In[14]:
+# In[7]:
 
 
 node.showControls()
 
 
-# In[17]:
+# In[21]:
 
 
 # Optional input, dictionary of Atom indices that should be matched in the search. 
 prematch = {}
 prematchstring = node.getInput("prematch")
-entries = prematchstring.split(",")
-for entry in entries:
-    idxA, idxB = entry.split("-")
-    prematch[ AtomIdx( int(idxA)) ] = AtomIdx( int(idxB) )
+if len(prematchstring) > 0: 
+    entries = prematchstring.split(",")
+    for entry in entries:
+        idxA, idxB = entry.split("-")
+        prematch[ AtomIdx( int(idxA)) ] = AtomIdx( int(idxB) )
+#print (prematch)
 
 
 # In[9]:
@@ -102,9 +116,12 @@ def writeLog(ligA, ligB, mapping):
 
 
 # Return a maximum of 10 matches, scored by RMSD and sorted from best to worst.
-mappings = BSS.Align.matchAtoms(lig1, lig2, matches=10, prematch=prematch)
+mappings = BSS.Align.matchAtoms(lig1, lig2, matches=10, prematch=prematch, scoring_function="RMSD", timeout=1*BSS.Units.Time.minute)
 # We retain the top mapping
 mapping = mappings[0]
+#print (len(mappings))
+#print (mappings)
+
 
 # In[13]:
 
@@ -120,7 +137,7 @@ def writeLog(ligA, ligB, mapping):
 system1.addMolecules(merged)
 
 
-# In[16]:
+# In[14]:
 
 
 # Log the mapping used
@@ -134,7 +151,7 @@ def writeLog(ligA, ligB, mapping):
 os.system(cmd)
 
 
-# In[20]:
+# In[15]:
 
 
 root = node.getInput("output")
@@ -145,22 +162,28 @@ def writeLog(ligA, ligB, mapping):
 mapping = "%s.mapping" % root
 
 
-# In[21]:
+# In[16]:
 
 
 cmd = "mv merged_at_lam0.pdb %s ; mv somd.pert %s ; mv somd.prm7 %s ; mv somd.rst7 %s ; mv somd.mapping %s ; rm somd.zip ; rm somd.cfg ; rm somd.err; rm somd.out" % (mergedpdb,pert,prm7,rst7,mapping)
 #print (cmd)
 os.system(cmd)
 
 
-# In[23]:
+# In[17]:
 
 
 node.setOutput("nodeoutput",[mergedpdb, pert, prm7, rst7, mapping])
 
 
-# In[24]:
+# In[18]:
 
 
 node.validate()
 
+
+# In[ ]:
+
+
+
+