add(obj.vina): prepare only once; prepare_each keyword is available t…

…o revert old behaviour

gaudi/algorithms.py

@@ -4,17 +4,17 @@
 ##############
 # GaudiMM: Genetic Algorithms with Unrestricted
 # Descriptors for Intuitive Molecular Modeling
-# 
+#
 # https://github.com/insilichem/gaudi
 #
 # Copyright 2017 Jaime Rodriguez-Guerra, Jean-Didier Marechal
-# 
+#
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
-# 
+#
 #      http://www.apache.org/licenses/LICENSE-2.0
-# 
+#
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -51,7 +51,7 @@
 
 
 def ea_mu_plus_lambda(population, toolbox, mu, lambda_, cxpb, mutpb, ngen, cfg,
-                      stats=None, halloffame=None, verbose=True, 
+                      stats=None, halloffame=None, verbose=True,
                       prompt_on_exception=True):
     """This is the :math:`(\mu + \lambda)` evolutionary algorithm.
 
@@ -112,7 +112,7 @@ def ea_mu_plus_lambda(population, toolbox, mu, lambda_, cxpb, mutpb, ngen, cfg,
     remaining_evals = estimated_evals - performed_evals
     remaining = timedelta(seconds=int(remaining_evals / speed))
     progress = '{:.2f}%'.format(100/(ngen+1))
-    logbook.record(gen=0, progress=progress, nevals=nevals, 
+    logbook.record(gen=0, progress=progress, nevals=nevals,
                    speed='{:.2f} ev/s'.format(speed), eta=remaining, **record)
     if verbose:
         logger.log(100, logbook.stream)
@@ -165,7 +165,7 @@ def ea_mu_plus_lambda(population, toolbox, mu, lambda_, cxpb, mutpb, ngen, cfg,
                     pass
             break
         else:
-            # Save a copy of an fully evaluated population, in case the 
+            # Save a copy of an fully evaluated population, in case the
             # simulation is stopped in next generation.
             population_ = population[:]
             if halloffame and cfg.output.check_every and not gen % cfg.output.check_every:

gaudi/base.py

@@ -4,17 +4,17 @@
 ##############
 # GaudiMM: Genetic Algorithms with Unrestricted
 # Descriptors for Intuitive Molecular Modeling
-# 
+#
 # https://github.com/insilichem/gaudi
 #
 # Copyright 2017 Jaime Rodriguez-Guerra, Jean-Didier Marechal
-# 
+#
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
-# 
+#
 #      http://www.apache.org/licenses/LICENSE-2.0
-# 
+#
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -100,7 +100,7 @@ def __init__(self, cfg=None, cache=None, **kwargs):
         self.expressed = False
         self.__ready__()
         self.__expression_hooks__()
-        
+
     def __ready__(self):
         """
         A *post-init* method used to avoid initialization problems with
@@ -237,7 +237,7 @@ def write(self, i, path=None):
         i : int
             Individual identificator in current generation or hall of fame
 
-        .. note :: 
+        .. note ::
 
             Maybe someday we can pickle it all :/
             >>> filename = os.path.join(path, '{}_{}.pickle.gz'.format(name,i))

gaudi/genes/torsion.py

@@ -4,17 +4,17 @@
 ##############
 # GaudiMM: Genetic Algorithms with Unrestricted
 # Descriptors for Intuitive Molecular Modeling
-# 
+#
 # https://github.com/insilichem/gaudi
 #
 # Copyright 2017 Jaime Rodriguez-Guerra, Jean-Didier Marechal
-# 
+#
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
-# 
+#
 #      http://www.apache.org/licenses/LICENSE-2.0
-# 
+#
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -43,7 +43,6 @@
 
 logger = logging.getLogger(__name__)
 
-
 def enable(**kwargs):
     kwargs = Torsion.validate(kwargs)
     return Torsion(**kwargs)
@@ -60,7 +59,7 @@ class Torsion(GeneProvider):
         Maximum number of degrees a bond can rotate
     max_bonds :
         Expected number of free rotations in molecule. Needed to store
-        arbitrary rotations. 
+        arbitrary rotations.
     anchor : str
         Molecule/atom_serial_number of reference atom for torsions
     rotatable_atom_types : list of str
@@ -69,7 +68,7 @@ class Torsion(GeneProvider):
     rotatable_atom_names : list of str
         Which type of atom names (as in chimera.Atom.name) should rotate.
         Defaults to ().
-    
+
     Attributes
     ----------
     allele : tuple of float
@@ -102,7 +101,7 @@ class Torsion(GeneProvider):
     BONDS_ROTS = {}
 
     def __init__(self, target=None, flexibility=360.0, max_bonds=None, anchor=None,
-                 rotatable_atom_types=('C3', 'N3', 'C2', 'N2', 'P'), 
+                 rotatable_atom_types=('C3', 'N3', 'C2', 'N2', 'P'),
                  rotatable_atom_names=(), rotatable_elements=(),
                  non_rotatable_bonds=(), precision=1, **kwargs):
         GeneProvider.__init__(self, **kwargs)
@@ -120,7 +119,7 @@ def __init__(self, target=None, flexibility=360.0, max_bonds=None, anchor=None,
 
     def __expression_hooks__(self):
         if self.max_bonds is None:
-            self.max_bonds = len(list(self.rotatable_bonds))
+            self.max_bonds = len(self.rotatable_bonds)
         self.allele = [self.random_angle() for i in xrange(self.max_bonds)]
 
     def express(self):
@@ -186,7 +185,7 @@ def rotatable_bonds(self):
 
         In this step, we have to discard non rotatable atoms (as requested
         by the user), and make sure the involved atoms are of compatible type.
-        Namely, one of them must be either C3, N3, C2 or N2, and both of them, 
+        Namely, one of them must be either C3, N3, C2 or N2, and both of them,
         non-terminal (more than one neighbor).
 
         If the bond is valid, get the BondRot object. Chimera will complain
@@ -204,7 +203,7 @@ def rotatable_bonds(self):
             return self.molecule._rotatable_bonds
 
     def _compute_rotatable_bonds(self):
-        bonds = sorted(self.molecule.bonds, 
+        bonds = sorted(self.molecule.bonds,
                        key=lambda b: min(y.serialNumber for y in b.atoms))
 
         non_rotatable_bonds = []
@@ -266,7 +265,7 @@ def anchor(self):
             else:
                 self.molecule._rotation_anchor = anchor
                 return anchor
-                
+
         target_gene = self.parent.find_molecule(self.target)
         try:
             search_gene = next(g for g in self.parent.genes.values()

gaudi/objectives/vina.py

@@ -56,6 +56,9 @@ class Vina(ObjectiveProvider):
         Key of the gene containing the molecule acting as receptor (protein)
     ligand : str
         Key of the gene containing the molecule acting as ligand
+    prepare_each : bool
+        Whether to prepare receptors and ligands in every evaluation or try
+        to cache the results for faster performance.
 
     Returns
     -------
@@ -65,13 +68,17 @@ class Vina(ObjectiveProvider):
     _validate = {
         parse.Required('receptor'): parse.Molecule_name,
         parse.Required('ligand'): parse.Molecule_name,
+        'prepare_each': bool,
         }
 
-    def __init__(self, receptor='Protein', ligand='Ligand', *args, **kwargs):
+    def __init__(self, receptor='Protein', ligand='Ligand', prepare_each=False,
+                 *args, **kwargs):
         ObjectiveProvider.__init__(self, **kwargs)
         self.receptor = receptor
         self.ligand = ligand
+        self.prepare_each = prepare_each
         self._paths = []
+        self._tmpfile = None
         if os.name == 'posix' and os.path.exists('/dev/shm'):
             self.tmpdir = '/dev/shm'
         else:
@@ -82,7 +89,6 @@ def evaluate(self, ind):
         Run a subprocess calling Vina binary with provided options,
         and parse the results. Clean tmp files at exit.
         """
-        self.tmpfile = os.path.join(self.tmpdir, next(_get_candidate_names()))
         receptor = ind.find_molecule(self.receptor)
         ligand = ind.find_molecule(self.ligand)
 
@@ -101,28 +107,57 @@ def evaluate(self, ind):
         finally:
             self.clean()
 
-    def prepare_receptor(self, molecule):
-        path = '{}_receptor.pdb'.format(self.tmpfile)
+    def _prepare(self, molecule, which='receptor'):
+        if which == 'receptor':
+            preparer = AD4ReceptorPreparation
+        elif which == 'ligand':
+            preparer = AD4LigandPreparation
+        else:
+            raise ValueError('which must be receptor or ligand')
+        path = '{}_{}.pdb'.format(self.tmpfile, which)
         pathqt = path + 'qt'
-        pdb = molecule.write(absolute=path, filetype='pdb')
-        self._paths.append(path)
-        mol = MolKit.Read(path)[0]
-        mol.buildBondsByDistance()
-        RPO = AD4ReceptorPreparation(mol, outputfilename=pathqt)
-        self._paths.append(pathqt)
+        if not os.path.isfile(pathqt):
+            pdb = molecule.write(absolute=path, filetype='pdb')
+            self._paths.append(path)
+            mol = MolKit.Read(path)[0]
+            mol.buildBondsByDistance()
+            RPO = preparer(mol, outputfilename=pathqt)
+            self._paths.append(pathqt)
+        else:
+            # update coordinates
+            self._update_pdbqt_coordinates(molecule.xyz(), pathqt)
         return pathqt
 
+    def prepare_receptor(self, molecule):
+        return self._prepare(molecule, 'receptor')
+
     def prepare_ligand(self, molecule):
-        path = '{}_ligand.pdb'.format(self.tmpfile)
-        pathqt = path + 'qt'
-        pdb = molecule.write(absolute=path, filetype='pdb')
-        self._paths.append(path)
-        mol = MolKit.Read(path)[0]
-        mol.buildBondsByDistance()
-        RPO = AD4LigandPreparation(mol, outputfilename=pathqt)
-        #    inactivate_all_torsions=True)
-        self._paths.append(pathqt)
-        return pathqt
+        return self._prepare(molecule, 'ligand')
+
+    @staticmethod
+    def _update_pdbqt_coordinates(xyz, path):
+        def is_atom(line):
+            if line[0:6] in ('ATOM  ', 'HETATM'):
+                for n in (line[30:38], line[38:46], line[46:54]):
+                    try:
+                        float(n)
+                    except:
+                        return False
+                return True
+            return False
+
+        with open(path, 'r+') as f:
+            lines = []
+            i = 0
+            for line in f:
+                if is_atom(line):
+                    line = line[:30] + '{:8.3f}{:8.3f}{:8.3f}'.format(*xyz[i]) + line[54:]
+                    i += 1
+                lines.append(line)
+            f.seek(0)
+            f.write(''.join(lines))
+            f.truncate()
+
 
     def parse_output(self, stream):
         for line in stream.splitlines():
@@ -131,6 +166,14 @@ def parse_output(self, stream):
         return -1000 * self.weight
 
     def clean(self):
+        if not self.prepare_each:
+            return
         for p in self._paths:
             os.remove(p)
-        self._paths = []
+        self._paths = []
+
+    @property
+    def tmpfile(self):
+        if self.prepare_each or self._tmpfile is None:
+            self._tmpfile = os.path.join(self.tmpdir, next(_get_candidate_names()))
+        return self._tmpfile