<reponame>insilichem/gpathfinder
#!/usr/bin/env python
# -*- coding: utf-8 -*-

##############
# GPathFinder: Identification of ligand pathways by a multi-objective
# genetic algorithm
#
# https://github.com/insilichem/gpathfinder
#
# Copyright 2019 <NAME>, <NAME>,
# <NAME>, <NAME>,
# <NAME> and <NAME>
#
# 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.
# See the License for the specific language governing permissions and
# limitations under the License.
##############

"""
GPathFinder allows to generate approximate ligand binding/unbinding pathways
to/from the binding site of a receptor (typically a protein).

This module must be used in conjunction with at least the following
three genes:
- A GPathFinder standard ``molecule`` gene for the ligand molecule.
- A GPathFinder standard ``molecule`` gene for the protein molecule.
- A GPathFinder ``path_torsion`` gene applied on the ligand to allow its
flexibility.

Optionally, the following genes can be used as well:
- A GPathFinder ``path_rotamers`` gene on the protein to allow conformational
changes in its side-chain.
- A GPathFinder``path_normal modes`` gene on the protein to allow the
exploration of its global folding through normal modes analysis.

To assess the quality of the pathways generated by this pool of genes,
the following objective must be used in the evaluation stage:
- At least one ``path scoring`` objective to evaluate how good are the
pathways generated by the genetic algorithm in terms of quality of the
sampling (avoid steric clashes, good vina score, good smoothness).
"""

# Python
from __future__ import print_function
import random
import logging
import numpy as np
from math import sqrt
import copy
from zipfile import ZipFile, ZIP_STORED
import os
import pprint

# Chimera
import chimera
import Matrix as M
from Measure import inertia
from Molecule import atom_positions
from chimera import Xform as X
from FitMap.search import random_rotation

# GPATH
from gpath.genes import GeneProvider
from gpath import parse
from gpath.exceptions import LigandDirectoryAndTorsion, ProteinDirectoryAndNM

IDENTITY = ((1.0, 0.0, 0.0, 0.0),
            (0.0, 1.0, 0.0, 0.0),
            (0.0, 0.0, 1.0, 0.0))

ZERO = chimera.Point(0.0, 0.0, 0.0)

pp = pprint.PrettyPrinter(4)
logger = logging.getLogger(__name__)

def enable(**kwargs):
    kwargs = Pathway.validate(kwargs)
    return Pathway(**kwargs)

class Pathway(GeneProvider):
    """
    Pathway class

    Parameters
    ----------
    ligand : str
        Name of the GPathFinder ``molecule`` gene containing the ligand
        molecule.
    protein : str
        Name of the GPathFinder ``molecule`` gene containing the receptor
        molecule.
    torsion_gene : str, optional
        Name of the GPathFinder ``path_torsion`` gene to allow flexibility of
        the ligand.
    rotamers_gene : str, optional
        Name of the GPathFinder ``path_rotamers`` gene to allow conformational
        changes on the protein side-chain.
    radius_rotamers: float, optional, defaults to 3.0
        Maximum distance from any point of the ligand in every frame
        that is searched for possible rotamers  of the protein side-chain.
    nm_gene : str, optional
        Name of the GPathFinder ``path_normalmodes`` gene to allow the
        exploration of the protein global folding through normal modes analysis.
    origin : 3-item list or tuple of float, optional
        Coordinates to a specific geometric center at which the origin
        of the ligand will be set, if they are different from the actual
        position of the ligand in the .mol2 file.
    inertia_axes : max of 6-item list or tuple of int, optional
        Starting point of the ligand when generating binding pathways
        will be set at one of the indicated ends of the inertia axes of
        the protein.
    destination : 3-item list or tuple of float, optional
        It indicates the coordinates to the geometric binding site or the
        expected end point of the pathway.
        If not set by the user, GPathFinder will assume an unbinding scenario,
        calculating automatically the distance from the ligand to outside
        the protein.
    max_step_separation : float, optional
        Maximum distance in Angstroms from one point of the ligand in
        the pathway to the next one. If not set by the user, GPathFinder
        calculates the value from the size of the ligand.
    min_step_increment : float, optional
        Minimum distance increment in Angstroms from the ligand's origin
        that has to be the ligand in one frame of the pathway with
        respect of the ligand's distance from the origin of the previous
        frame of the pathway. If not set by the user, GPathFinder calculates
        the value as 1/5 by max_step_separation.
    frontier_margin : float, optional, defaults to 0.0
        Safety distance from the frontier of the protein to consider that
        the ligand is outside. Specially useful when using large ligands
        that can be stucked at the frontier.
    mut_pos_pb : float, optional, defaults to 0.10
        When a mutation occurs, this value is the probability of such
        mutation to be of the type `positions`, that is, the mutation
        changes the actual trajectory of the pathway. Warning: parameter
        for advanced users, usually the default value is correct for the
        vast majority of the systems.

    Attributes
    ----------
    allele : dict
        allele['positions'] = list of translation matrices of the ligand
                              along the frames of the pathway, as
                              explained in Notes.
        allele['rotations'] = list of rotation materices of the ligand
                              along the frames of the pathway, as
                              explained in Notes.
        allele['torsions'] = list of torsion alleles applied to the
                             ligand along the frames of the pathway (see
                             allele of GPathFinder ``path_torsion`` gene).
        allele['torsion_anchor'] = [str, int] ligand gene name and
                                   serial number of the anchor atom used
                                   in the torsion gene.
        allele['rotamers'] = list of data of the rotamers applied on the
                             protein along the frames of the pathway,
                             as explained in Notes.
        allele['normal_modes'] = list of normal modes alleles applied to
                                 the protein along the frames of the
                                 pathway (see allele of GPathFinder
                                 ``path_normalmodes`` gene).
    scores : list of dict
        Scores reported for every frame of the pathway. Each key of the
        dict correspond to one score (e.g. clashes).

    Notes
    -----
    ** Priority and configuration of the origin point of the ligand **

    The origin point/s of the ligand can be provided by three different
    manners:

        1. The ends of the inertia axes of the protein provided with the
           ``inertia_axes``.
        2. The coordinates set in ``origin`` parameter.
        3. The actual ubication of the ligand in its .mol2 file.

    The priority in case of conflict will be the order listed before.
    That is, if ``ìnertia_axes`` parameter is provided, this will
    override a possible ``origin`` parameter and the actual ubication in
    the .mol2 file.

    ** Details of the allele structure **

    - ``positions`` of the ligand along the pathway will be defined by
      translation matrices with the following shape:

        (
        (1, 0, 0, Tx),
        (0, 1, 0, Ty),
        (0, 0, 1, Tz)
        )

    - ``rotations`` of the ligand along the pathway will be defined by
      rotation matrices with the following shape:

        (
        (R1, R2, R3, 0),
        (R4, R5, R6, 0),
        (R7, R8, R9, 0)
        )

    - ``rotamers`` applied to the ligand at one frame of the pathway
      will be defined as 3 lists, each one containing:

        1. Residue numbers over which the rotamers are applied.
        2. Actual rotamer allele (see allele of GPathFinder ``path_rotamers`` gene)
        3. Chi angles corresponding to each rotamer.
    """

    _validate = {
        parse.Required('ligand'): parse.Molecule_name,
        parse.Required('protein'): parse.Molecule_name,
        'torsion_gene': parse.Molecule_name,
        'rotamers_gene': parse.Molecule_name,
        'radius_rotamers': parse.All(parse.Coerce(float), parse.Range(min=0)),
        'nm_gene': parse.Molecule_name,
        'origin': parse.Coordinates,
        'inertia_axes': [parse.All(parse.Coerce(int),
                            parse.Range(min=0, max=5))],
        'destination': parse.Coordinates,
        'max_step_separation': parse.Coerce(float),
        'min_step_increment': parse.Coerce(float),
        'frontier_margin': parse.Coerce(float),
        'mut_pos_pb': parse.Coerce(float),
        }

    def __init__(self, ligand, protein, torsion_gene=None, rotamers_gene=None,
                radius_rotamers=3.0, nm_gene=None,  origin=None,
                inertia_axes=None, destination=None, max_step_separation=None,
                min_step_increment=None, frontier_margin=0.0, mut_pos_pb=0.1, **kwargs):
        GeneProvider.__init__(self, **kwargs)
        self.ligand = ligand
        self.protein = protein
        self.torsion_gene = torsion_gene
        self.rotamers_gene = rotamers_gene
        self.radius_rotamers = radius_rotamers
        self.nm_gene = nm_gene
        self.destination = destination
        self.origin = origin if origin else None
        self.inertia_axes = inertia_axes if inertia_axes else None
        self.max_separation = max_step_separation
        self.min_increment = min_step_increment
        self.frontier_margin = frontier_margin
        self.mut_pos_pb = mut_pos_pb
        self.act_rotamers = []
        self.which_evaluations = []

    def __ready__(self):
        # Check if ligand conformers and torsion are used simultaneously
        if len(self.ligand_g.catalog) > 1 and self.torsion_gene:
            raise LigandDirectoryAndTorsion('If you set the Ligand molecule to be a directory, you can not use a path_torsion gene simultaneously')
        # Check if protein conformers and NM are used simultaneously
        if len(self.protein_g.catalog) > 1 and self.nm_gene:
            raise ProteinDirectoryAndNM('If you set the Protein molecule to be a directory, you can not use a path_normalmodes gene simultaneously')

        # Calculate inertia axes of the protein for further use
        self.p_axes, self.p_d2, self.p_center = inertia.atoms_inertia(self.protein_mol.atoms)
        self.p_elen = [a + self.frontier_margin for a in inertia.inertia_ellipsoid_size(self.p_d2)]
        self.p_axes = np.array(self.p_axes)

        if self.inertia_axes: #Prepare origin from inertia axes
            origin_points = []
            for axis, length in zip(self.p_axes, self.p_elen):
                origin_points.append((axis*length)+self.p_center)
                origin_points.append((-axis*length)+self.p_center)
            points = []
            for i, point in enumerate(origin_points):
                if i in self.inertia_axes:
                    points.append(point)
            self.origin = random.choice(points)
        elif not self.origin:
            #Origin is the geometric center of the ligand if not
            #provided explicitly
            self.origin = self.ligand_center

        #Calculate min_increment and max_separation if not set by user
        if not self.max_separation:
            axes, d2, center = inertia.atoms_inertia(self.ligand_mol.atoms)
            elen = [a for a in inertia.inertia_ellipsoid_size(d2)]
            self.max_separation = min(elen) if min(elen) > 1.0 else 1.0
        if self.min_increment is None:
            self.min_increment = 1.0 * self.max_separation / 5.0

        #Operations permitted to improve the path in mutate
        self._op_roulette = ['rotation']
        if self.torsion_gene:
            self._op_roulette.append('torsion')
        if self.rotamers_gene:
            self._op_roulette.append('rotamers')
        if self.nm_gene:
            self._op_roulette.append('nm')
        if len(self.protein_g.catalog) > 1:
            self._op_roulette.append('protein') #Conformers for the protein
        if len(self.ligand_g.catalog) > 1:
            self._op_roulette.append('ligand') #Conformers for the ligand

        #Making the initial allele
        self.allele = {}
        self.allele['positions'] = [((1.0, 0.0, 0.0, self.origin[0]),
                                    (0.0, 1.0, 0.0, self.origin[1]),
                                    (0.0, 0.0, 1.0, self.origin[2]))]
        self.allele['rotations'] = [IDENTITY]
        self.allele['min_increment'] = self.min_increment
        self.allele['max_separation'] = self.max_separation
        if self.torsion_gene:
            self.allele['torsions'] = [[]]
            self.allele['torsion_anchor'] = list(self.torsion_g._anchor)
        if self.rotamers_gene:
            self.allele['rotamers'] = [ [[],[],[]] ]
        if self.nm_gene:
            self.allele['normal_modes'] = [None]
        self.allele['protein'] = [0] #Conformer for the protein
        self.allele['ligand'] = [0] #Conformer for the ligand
        self.allele['coord_residues'] = [[]] #To store possible coord residues when evaluating metal_sites
        self.allele['mate_torsions'] = 0  #To register if mate is being useful
        self.allele['mate_positions'] = 0 #To register if mate is being useful
        self.allele['mate_rotations'] = 0 #To register if mate is being useful
        self.allele['mutate_positions'] = 0 #To register if mutate is being useful
        self.allele['mutate_torsions'] = 0 #To register if mutate is being useful
        self.allele['mutate_rotations'] = 0 #To register if mutate is being useful
        self.allele['mutate_rotamers'] = 0 #To register if mutate is being useful
        self.allele['mutate_nm'] = 0 #To register if mutate is being useful
        self.allele['mutate_protein'] = 0 #To register if mutate is being useful
        self.allele['mutate_ligand'] = 0 #To register if mutate is being useful
        self.scores = [{}] #To store scores at the evaluation stage

        #Initial pathway
        last_point = [x[3] for x in self.allele['positions'][-1]]
        if self.destination:
            not_arrived = distance(self.destination, last_point) > (self.max_separation)
        else:
            y = np.array([a-b for a,b in zip(last_point, self.p_center)])
            y = np.absolute(self.p_axes.dot(y))
            not_arrived = is_inside_ellipsoid(self.p_elen, y)

        while not_arrived:
            new_point = random_extension_point(last_point, self.max_separation,
                                self.origin, min_increment=self.min_increment,
                                destination=self.destination)
            if new_point is not None:
                self.allele['positions'].append(new_point)
                self.allele['rotations'].append(random_rotation())
                if self.torsion_gene:
                    self.torsion_g.gp_mutate(1.0)
                    self.allele['torsions'].append(copy.deepcopy(self.torsion_g.allele))
                if self.nm_gene:
                    self.allele['normal_modes'].append(random.randint(0,
                                    self.parent.genes[self.nm_gene].n_samples))
                self.allele['protein'].append(random.randint(0, len(self.protein_g.catalog)-1))
                self.allele['ligand'].append(random.randint(0, len(self.ligand_g.catalog)-1))
                if self.rotamers_gene:
                    self.allele['rotamers'].append([[],[],[]])
                self.allele['coord_residues'].append([])
                self.scores.append({})
                last_point = [x[3] for x in self.allele['positions'][-1]]
            if self.destination:
                not_arrived = distance(self.destination, last_point) > (self.max_separation)
            else:
                y = np.array([a-b for a,b in zip(last_point, self.p_center)])
                y = np.absolute(self.p_axes.dot(y))
                not_arrived = is_inside_ellipsoid(self.p_elen, y)

    def express(self):
        """
        Expression is controlled in gp_express to avoid unnecessary
        calls and supervise what frames are expressed.
        """
        pass

    def unexpress(self):
        """
        Unexpression is controlled in gp_unexpress to avoid unnecessary
        calls and supervise what frames are expressed.
        """
        pass

    def gp_express(self, i, with_rotamers=True, smoothness=False):
        """
        For the demanded frame ``i``:
            1. Apply translation, rotation and torsions to the ligand.
            2. Apply, if present, rotamers and normal modes to the
               protein.
        Optional parameter ``ẁith_rotamers`` controls if rotamers have
        to be expressed or not (useful when a rotamer actualization is
        pending).
        Optional parameter ``smoothness`` is used to express only the
        ligand transformations in order to calculate RMSD between frames.
        If True, only torsions and rotations of the ligand are applied.
        """
        ligand_sample = self.allele['ligand'][i]
        self.ligand_g.allele = self.ligand_g.catalog[ligand_sample]
        self.ligand_g._need_express = True
        self.ligand_g.express()

        protein_sample = self.allele['protein'][i]
        self.protein_g.allele = self.protein_g.catalog[protein_sample]
        self.protein_g._need_express = True
        self.protein_g.express()

        or_x, or_y, or_z = self.ligand_center
        to_zero = ((1.0, 0.0, 0.0, -or_x),
                    (0.0, 1.0, 0.0, -or_y),
                    (0.0, 0.0, 1.0, -or_z))

        rotations = self.allele['rotations']
        if smoothness:
            matrices = (IDENTITY,) + (rotations[i],) + (to_zero,)
        else:
            positions = self.allele['positions']
            matrices = (positions[i],) + (rotations[i],) + (to_zero,)
        matrices = M.multiply_matrices(*matrices)
        self.ligand_mol.openState.xform = M.chimera_xform(matrices)

        if self.torsion_gene:
            self.torsion_g.allele = copy.deepcopy(self.allele['torsions'][i])
            self.torsion_g._need_express = True
            self.torsion_g.express()

        if self.nm_gene and not smoothness:
            sample_number = self.allele['normal_modes'][i]
            if sample_number:
                self.nm_g.allele = self.nm_g.NORMAL_MODES_SAMPLES[sample_number]
                self.nm_g._need_express = True
                self.nm_g.express()

        if self.rotamers_gene and with_rotamers and not smoothness:
            if i in self.act_rotamers:
                #Actualize rotamers
                self.actualize_rotamers(i)

            if self.allele['rotamers'][i][1]:
                residues = []
                for rot in self.allele['rotamers'][i][0]:
                    a, b = rot.split('/')
                    residues.append([a, int(b)])

                self.rotamers_g._residues = residues
                self.rotamers_g.__ready__()
                self.rotamers_g.allele = self.allele['rotamers'][i][1]
                for (molname, pos), residue in self.rotamers_g.residues.items():
                    self.rotamers_g.allele[residues.index([molname, pos])] = self.allele['rotamers'][i][1][residues.index([molname, pos])]
                self.rotamers_g._need_express = True
                self.rotamers_g.express()

    def gp_unexpress(self, i, smoothness=False):
        """
        For the demanded frame ``i``:
            1. Undo translation, rotation and torsions made to the
               ligand.
            2. Undo, if present, rotamers and normal modes changes made
               to the protein.
        Optional parameter ``smoothness`` is used to unexpress only the
        ligand trasnformations when used to calculate RMSD of the ligand
        between frames. If True, only torsions and rotations of the
        ligand are unapplied.
        """
        self.ligand_mol.openState.xform = X()

        if self.torsion_gene:
            self.torsion_g.unexpress()
            self.torsion_g._need_express = False

        if self.rotamers_gene and not smoothness:
            if self.allele['rotamers'][i][0]:
                self.rotamers_g.unexpress()
                self.rotamers_g._need_express = False

        if self.nm_gene and not smoothness:
            sample_number = self.allele['normal_modes'][i]
            if sample_number:
                self.nm_g.unexpress()
                self.nm_g._need_express = False

        self.protein_g.unexpress()
        self.protein_g._need_express = False

        self.ligand_g.unexpress()
        self.ligand_g._need_express = False

    def mate(self, mate):
        """
        Randomnly choose an operation between [translations, rotations,
        torsions] of the ligand and make a mate of the allele of the two
        individuals.
        """
        #frame for self is the selected by score probability
        k = random.choice(self.scores[0].keys())
        scores_list = [sc[k] for sc in self.scores]
        m = min(scores_list)
        if m < 0:
            scores_list = [sc - m for sc in scores_list]
        scores_list = [sc + 1.0 for sc in scores_list] #to avoid divisions by 0 when calculating probs
        probs = [float(sc)/sum(scores_list[1:]) for sc in scores_list[1:]]
        i = np.random.choice(range(1, len(self.allele['positions'])), p=probs)

        #frame for mate is the nearest to the selected for self
        self_pos = [x[3] for x in self.allele['positions'][i]]
        dists = []
        for frame in range(1, len(mate.allele['positions'])):
            mate_pos = [x[3] for x in mate.allele['positions'][frame]]
            dists.append(distance(self_pos, mate_pos))
        j = dists.index(min(dists)) + 1 #dists doesn't contain distance for the frame 0

        if self.torsion_gene:
            op = random.choice([self.mate_torsions, self.mate_rotations])
        else:
            op = self.mate_rotations

        op(mate, i, j)

    def mate_rotations(self, mate, i, j):
        """
        For the frame ``i`` of self individual and ``j`` of the mate
        individual, make a interpolation of the rotations of the two
        ligands and assign the resulting matrix to both.
        """
        self.allele['mate_rotations'] = self.allele['mate_rotations'] + 1  #For register if mate is being useful
        mate.allele['mate_rotations'] = mate.allele['mate_rotations'] + 1  #For register if mate is being useful

        interp = interpolate_rotations(self.allele['rotations'][i],
                                            mate.allele['rotations'][j])
        self.allele['rotations'][i] = copy.deepcopy(interp)
        mate.allele['rotations'][j] = copy.deepcopy(interp)
        #Delete previous calculated scores
        self.scores[i] = {}
        mate.scores[j] = {}
        #Force rotamers actualization
        self.act_rotamers.append(i)
        mate.act_rotamers.append(j)

    def mate_torsions(self, mate, i, j):
        """
        For the frame ``i`` of self individual and ``j`` of the mate
        individual, make a torsion gp_mate (see ``path_torsion`` gene).
        """
        self.allele['mate_torsions'] = self.allele['mate_torsions'] + 1  #For register if mate is being useful
        mate.allele['mate_torsions'] = mate.allele['mate_torsions'] + 1  #For register if mate is being useful

        self.torsion_g.allele = copy.deepcopy(self.allele['torsions'][i])
        mate.torsion_g.allele = copy.deepcopy(mate.allele['torsions'][j])
        self.torsion_g.gp_mate(mate.torsion_g)
        self.allele['torsions'][i] = copy.deepcopy(self.torsion_g.allele)
        mate.allele['torsions'][j] = copy.deepcopy(mate.torsion_g.allele)
        #Delete previous calculated scores
        self.scores[i] = {}
        mate.scores[j] = {}
        #Force rotamers actualization
        self.act_rotamers.append(i)
        mate.act_rotamers.append(j)

    def mutate(self, indpb):
        """
        There are two possible mutations:
            1. Cut the pathway and remake  all the frames from the
               selected cutting frame (with probability mut_pos_pb).
            2. Change the [rotation, torsions, rotamers, normal modes]
               of one frame.
        The frame affected by the mutation is selected randomnly, with
        the probabilities of each frame weighted by their scores (more
        probability to be mutated if the frame has a worse score).
        The operation between the available pool is chosen randomly .
        """
        if random.random() < self.indpb:
            k = random.choice(self.scores[0].keys())
            scores_list = [sc[k] for sc in self.scores]
            m = min(scores_list)
            if m < 0:
                scores_list = [sc - m for sc in scores_list]
            scores_list = [sc + 1.0 for sc in scores_list] #to avoid divisions by 0 when calculating probs
            probs = [float(sc)/sum(scores_list[1:]) for sc in scores_list[1:]]
            i = np.random.choice(range(1, len(self.allele['positions'])), p=probs)

            n = random.random()
            if n < self.mut_pos_pb: #Modify the positions of the pathway
                i = np.random.choice(range(1, len(self.allele['positions']))) # No prob biass in this case
                self.allele['mutate_positions'] = self.allele['mutate_positions'] + 1 #For register if mutate is being useful
                #Delete frames from i
                self.allele['positions'] = self.allele['positions'][:i]
                self.allele['rotations'] = self.allele['rotations'][:i]
                if self.torsion_gene:
                    self.allele['torsions'] = self.allele['torsions'][:i]
                if self.rotamers_gene:
                    self.allele['rotamers'] = self.allele['rotamers'][:i]
                if self.nm_gene:
                    self.allele['normal_modes'] = self.allele['normal_modes'][:i]
                self.allele['protein'] = self.allele['protein'][:i]
                self.allele['ligand'] = self.allele['ligand'][:i]
                self.allele['coord_residues'] = self.allele['coord_residues'][:i]
                self.scores = self.scores[:i]
                #Generate new frames
                last_point = [x[3] for x in self.allele['positions'][-1]]
                if self.destination:
                    not_arrived = distance(self.destination, last_point) > (self.max_separation)
                else:
                    y = np.array([a-b for a,b in zip(last_point, self.p_center)])
                    y = np.absolute(self.p_axes.dot(y))
                    not_arrived = is_inside_ellipsoid(self.p_elen, y)

                while not_arrived:
                    new_point = random_extension_point(last_point, self.max_separation,
                                        self.origin, min_increment=self.min_increment,
                                        destination=self.destination)
                    if new_point is not None:
                        self.allele['positions'].append(new_point)
                        self.allele['rotations'].append(random_rotation())
                        if self.torsion_gene:
                            self.torsion_g.gp_mutate(1.0)
                            self.allele['torsions'].append(copy.deepcopy(self.torsion_g.allele))
                        if self.nm_gene:
                            self.allele['normal_modes'].append(random.randint(0,
                                            self.parent.genes[self.nm_gene].n_samples))
                        self.allele['protein'].append(random.randint(0, len(self.protein_g.catalog)-1))
                        self.allele['ligand'].append(random.randint(0, len(self.ligand_g.catalog)-1))
                        if self.rotamers_gene:
                            self.allele['rotamers'].append([[],[],[]])
                        self.allele['coord_residues'].append([])
                        self.scores.append({})
                        last_point = [x[3] for x in self.allele['positions'][-1]]
                    if self.destination:
                        not_arrived = distance(self.destination, last_point) > (self.max_separation)
                    else:
                        y = np.array([a-b for a,b in zip(last_point, self.p_center)])
                        y = np.absolute(self.p_axes.dot(y))
                        not_arrived = is_inside_ellipsoid(self.p_elen, y)
            else:
                #Try to modify rotation/torsion/rotamers/nm/sample protein
                operation = random.choice(self._op_roulette)
                if operation == 'rotation':
                    self.allele['mutate_rotations'] = self.allele['mutate_rotations'] + 1 #For register if mutate is being useful
                    self.allele['rotations'][i] = random_rotation()
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'torsion':
                    self.allele['mutate_torsions'] = self.allele['mutate_torsions'] + 1 #For register if mutate is being useful
                    self.torsion_g.gp_mutate(1.0)
                    self.allele['torsions'][i] = copy.deepcopy(self.torsion_g.allele)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'rotamers':
                    self.allele['mutate_rotamers'] = self.allele['mutate_rotamers'] + 1 #For register if mutate is being useful
                    self.scores[i] = {} #put the score to 0 and allow path_scoring to recalculate rotamers
                     #Search which rotamers are inside the search radius
                    residues = set()
                    atoms = [a for a in surrounding_atoms(self.ligand_mol, self.protein_mol, self.radius_rotamers)]
                    for atom in atoms:
                        residues.add(atom.residue.id.position)
                    residues_list = ['{}/{}'.format(self.protein, r) for r in residues]
                    residues = []
                    for rot in residues_list:
                        a, b = rot.split('/')
                        residues.append([a, int(b)])
                    #Create the new rotamers allele and express it to store in path allele
                    self.rotamers_g._residues = residues
                    self.rotamers_g.__ready__()
                    self.rotamers_g._need_express = True
                    self.rotamers_g.express()
                    #Store the new rotamer data in pathway gene
                    self.allele['rotamers'][i][0] = []
                    self.allele['rotamers'][i][1] = []
                    self.allele['rotamers'][i][2] = []
                    for ((molname, pos), residue), a in \
                                                    zip(self.rotamers_g.residues.items(),
                                                    self.rotamers_g.allele):
                        res_name = str(molname)+'/'+str(pos)
                        self.allele['rotamers'][i][0].append(res_name)
                        self.allele['rotamers'][i][1].append(a)
                        self.allele['rotamers'][i][2].append(self.rotamers_g.all_chis(residue))
                elif operation == 'nm':
                    self.allele['mutate_nm'] = self.allele['mutate_nm'] + 1 #For register if mutate is being useful
                    self.allele['normal_modes'][i] = random.randint(0,
                                    self.parent.genes[self.nm_gene].n_samples)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'protein':
                    self.allele['mutate_protein'] = self.allele['mutate_protein'] + 1 #For register if mutate is being useful
                    self.allele['protein'][i] = random.randint(0, len(self.protein_g.catalog)-1)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'ligand':
                    self.allele['mutate_ligand'] = self.allele['mutate_ligand'] + 1 #For register if mutate is being useful
                    self.allele['ligand'][i] = random.randint(0, len(self.ligand_g.catalog)-1)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)

    def actualize_rotamers(self, i):
        """
        Actualize the rotamers list and allele when a position/
        orientation/torsions of a ligand have changed (it could be
        possible that with the new ligand expression the surrounding
        rotamers would be different, so it is necessary to actualize the
        information). Old rotamer allele is maintained when possible
        (i.e. the rotamers that are present in the previous and the new
        sourrondings of the ligand).
        """
        #Search which rotamers are inside the search radius
        residues = set()
        atoms = [a for a in surrounding_atoms(self.ligand_mol, self.protein_mol, self.radius_rotamers)]
        for atom in atoms:
            residues.add(atom.residue.id.position)
        residues_list = ['{}/{}'.format(self.protein, r) for r in residues]
        residues = []
        for rot in residues_list:
            a, b = rot.split('/')
            residues.append([a, int(b)])

        #Create the new rotamers allele and express i
        self.rotamers_g._residues = residues
        self.rotamers_g.__ready__()
        self.rotamers_g._need_express = True
        self.rotamers_g.express()

        #Store the new rotamer data in pathway gene (only new residues)
        old_res = [[],[],[]]
        old_res[0] = list(self.allele['rotamers'][i][0])
        old_res[1] = list(self.allele['rotamers'][i][1])
        old_res[2] = list(self.allele['rotamers'][i][2])
        self.allele['rotamers'][i][0] = []
        self.allele['rotamers'][i][1] = []
        self.allele['rotamers'][i][2] = []
        new_res = [str(molname)+'/'+str(pos) for ((molname, pos), residue), a in \
                    zip(self.rotamers_g.residues.items(),
                        self.rotamers_g.allele)]
        #Incorporate old residues present in the new region as is
        for j, res_name in enumerate(old_res[0]):
            if res_name in new_res:
                self.allele['rotamers'][i][0].append(res_name)
                self.allele['rotamers'][i][1].append(old_res[1][j])
                self.allele['rotamers'][i][2].append(old_res[2][j])
        #Incorporate new residues of the new region
        for ((molname, pos), residue), a in \
                    zip(self.rotamers_g.residues.items(),
                        self.rotamers_g.allele):
            res_name = str(molname)+'/'+str(pos)
            if res_name not in self.allele['rotamers'][i][0]:
                self.allele['rotamers'][i][0].append(res_name)
                self.allele['rotamers'][i][1].append(a)
                self.allele['rotamers'][i][2].append(self.rotamers_g.all_chis(residue))

        self.rotamers_g.unexpress()
        self.rotamers_g._need_express = False

        self.act_rotamers.remove(i) #Rotamers are actualized

    def write(self, path, name, *args, **kwargs):
        """
        Each individual output is a .zip file that contains:
            1. A .txt file with the allele data (allele.txt).
            2. A .txt file with the scoring data (scores.txt).
            3. One .pdb file for each frame of the calculated pathway
               (named frame_xxx.pdb, where xxx is the number of frame).
        """
        fullname = os.path.join(path, '{}_{}.zip'.format(name, self.name))
        fullname = self.create_output_zip(path, fullname)
        return fullname

    def create_output_zip(self, path, fullname):
        #Txt file with the allele
        allelename = os.path.join(path, 'allele.txt')
        with open(allelename, 'w') as f:
            f.write(pp.pformat(self.allele))

        #Txt file with the scores
        scoresname = os.path.join(path, 'scores.txt')
        with open(scoresname, 'w') as f:
            f.write(pp.pformat(self.scores))

        #Pdb file with the trajectory of the pathway
        trajname = os.path.join(path, 'trajectory.pdb')
        points = self.allele['positions']
        mol = chimera.Molecule()
        r = mol.newResidue("path", " ", 1, " ")
        atoms = []
        atom = mol.newAtom("0", chimera.Element("H"))
        atoms.append(atom)
        point = [x[3] for x in points[0]]
        atoms[-1].setCoord(chimera.Point(*point))
        r.addAtom(atoms[-1])
        for i in range(1, len(points)):
            point = [x[3] for x in points[i]]
            atom = mol.newAtom(str(i), chimera.Element("H"))
            atoms.append(atom)
            atoms[-1].setCoord(chimera.Point(*point))
            r.addAtom(atoms[-1])
            mol.newBond(atoms[i-1], atoms[i])
        chimera.pdbWrite([mol], chimera.Xform(), trajname)

        with ZipFile(fullname, 'w', ZIP_STORED) as z:
            z.write(allelename, os.path.basename(allelename))
            os.remove(allelename)
            z.write(scoresname, os.path.basename(scoresname))
            os.remove(scoresname)
            z.write(trajname, os.path.basename(trajname))
            os.remove(trajname)

            #Pdb files with the actual frames of the pathway
            for i in range(0, len(self.allele['positions'])):

                self.gp_express(i)
                framename = os.path.join(path, "frame_{:03d}.pdb".format(i))
                chimera.pdbWrite([self.protein_mol, self.ligand_mol],
                                    chimera.Xform(), framename)

                # Prepare pdb files with format required by Chimera MD movie import
                with open(framename) as frame_file:
                    pdb_file_lines = []
                    line = frame_file.readline()
                    while line:
                        if line[:4] == "ATOM" or line[:6] == "HETATM":
                            pdb_file_lines.append(line)
                        line = frame_file.readline()

                os.remove(framename)
                with open(framename, 'w') as f:
                    prev_atom = 0
                    prev_res = 0
                    last_change_res = 0
                    for line in pdb_file_lines:
                        current_atom, current_res = int(line[6:11]), int(line[22:26])
                        if current_atom > prev_atom:
                            prev_atom = current_atom
                        else:
                            current_atom = prev_atom + 1
                            prev_atom = current_atom
                        if current_res > prev_res:
                            prev_res = current_res
                            last_change_res = current_res
                        elif current_res != last_change_res:
                            last_change_res = current_res
                            prev_res += 1
                            current_res = prev_res
                        else:
                            current_res = prev_res
                        line = line[:6] + '{:5d}'.format(current_atom) + line[11:22] + '{:4d}'.format(current_res) + line[26:]
                        f.write("%s" % line)

                z.write(framename, os.path.basename(framename))
                os.remove(framename)
                self.gp_unexpress(i)

        return fullname

    #####
    @property
    def ligand_mol(self):
        return self.parent.find_molecule(self.ligand).compound.mol

    @property
    def protein_mol(self):
        return self.parent.find_molecule(self.protein).compound.mol

    @property
    def ligand_g(self):
        return self.parent.find_molecule(self.ligand)

    @property
    def protein_g(self):
        return self.parent.find_molecule(self.protein)

    @property
    def torsion_g(self):
        if self.torsion_gene:
            return self.parent.genes[self.torsion_gene]
        else:
            return None

    @property
    def rotamers_g(self):
        if self.rotamers_gene:
            return self.parent.genes[self.rotamers_gene]
        else:
            return None

    @property
    def nm_g(self):
        if self.nm_gene:
            return self.parent.genes[self.nm_gene]
        else:
            return None

    @property
    #Geometric center of the ligand
    def ligand_center(self):
        coordinates = atom_positions(self.ligand_mol.atoms)
        c = np.average(coordinates, axis=0)
        return c


#############
# Some useful functions
def random_extension_point(center, r, origin, min_increment=0, destination=None):
    """
    Calculate a new random coordinates that acomplishes the following
    conditions:
        1. The distance from the ``center`` has to be <= r.
        2. The distance from the origin has to be >= than that from the
           ``center`` to origin plus ``min_increment`` when calculating
           pathways without a ``destination``.
        3. If a ``destination`` is set, the new point has to be nearer
           to this destination than the ``center`` point.

    Parameters
    ----------
    center : 3-tuple of float
        Coordinates of the previous point (origin of the search sphere).
    r : float
        Radius of the search sphere (maximum distance at will the new
        point will be from ``center``).
    origin : 3-tuple of float
        Coordinates of the origin point of the pathway.
    min_increment : float, optional, defaults to 0
        Minimum increment of distance from the origin / to the destination
        of the new point respect to th previous one.
    destination : float or 3-tuple of float, optional
        Distance from the origin or coordinates of the destination point
        of the pathway.

    Returns
    -------
    A 3-tuple of float representing the coordinates that accomplish all
    the criteria. If a valid point is not obtained in 100 iterations,
    returns None.

    Notes
    -----
    The combination of the ``r`` (max_step_separation) and
    ``min_increment`` will control the distance of the different frames
    of the pathway between them and the possibility of deviation from the
    straight line (as higher the difference between these parameters,
    more deviation from the straight line will be allowed).
    """
    for i in range(0, 100):
        new_point = [random.uniform(-r, r) for m in center]
        new_point = [x + y for x, y in zip(new_point, center)]
        x, y, z = new_point
        if not isinstance(destination, list):
            if ((distance(new_point, origin) >= (distance(center, origin) +
                                                min_increment))
                and (distance(new_point, center) <= r)):
                return ((1.0, 0.0, 0.0, x),
                        (0.0, 1.0, 0.0, y),
                        (0.0, 0.0, 1.0, z))
        else:
            if (((distance(new_point, destination) +
                    min_increment) < distance(center, destination))
                and (distance(new_point, center) <= r)):
                return ((1.0, 0.0, 0.0, x),
                        (0.0, 1.0, 0.0, y),
                        (0.0, 0.0, 1.0, z))
    return None

def distance(a, b):
    ax, ay, az = a
    bx, by, bz = b
    return(sqrt((bx-ax)*(bx-ax)+(by-ay)*(by-ay)+(bz-az)*(bz-az)))

def interpolate_rotations(a, b):
    xf1 = M.chimera_xform(a)
    xf2 = M.chimera_xform(b)
    interp = M.xform_matrix(M.interpolate_xforms(xf1, ZERO, xf2, 0.5))
    return interp

def is_inside_ellipsoid(elen, point):
    discr = (point[0]/elen[0])**2 + (point[1]/elen[1])**2 + (point[2]/elen[2])**2
    if discr < 1.0:
        return True
    else:
        return False

def surrounding_atoms(ligand, protein, radius_rotamers):
    """
    Get atoms in the search zone, based on the molecule and the
    radius.
    """
    z = chimera.selection.ItemizedSelection()
    z.add([a for a in ligand.atoms])
    z.merge(chimera.selection.REPLACE,
                    chimera.specifier.zone(z, 'atom', None,
                                            radius_rotamers,
                                            [ligand, protein]))
    atom_list = [at for at in z.atoms() if at not in ligand.atoms]
    return atom_list

#!/usr/bin/env python
# -*- coding: utf-8 -*-

##############
# GPathFinder: Identification of ligand pathways by a multi-objective
# genetic algorithm
#
# https://github.com/insilichem/gpathfinder
#
# Copyright 2019 <NAME>, <NAME>,
# <NAME>, <NAME>,
# <NAME> and <NAME>
#
# 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.
# See the License for the specific language governing permissions and
# limitations under the License.
##############

"""
GPathFinder allows to generate approximate ligand binding/unbinding pathways
to/from the binding site of a receptor (typically a protein).

This module must be used in conjunction with at least the following
three genes:
- A GPathFinder standard ``molecule`` gene for the ligand molecule.
- A GPathFinder standard ``molecule`` gene for the protein molecule.
- A GPathFinder ``path_torsion`` gene applied on the ligand to allow its
flexibility.

Optionally, the following genes can be used as well:
- A GPathFinder ``path_rotamers`` gene on the protein to allow conformational
changes in its side-chain.
- A GPathFinder``path_normal modes`` gene on the protein to allow the
exploration of its global folding through normal modes analysis.

To assess the quality of the pathways generated by this pool of genes,
the following objective must be used in the evaluation stage:
- At least one ``path scoring`` objective to evaluate how good are the
pathways generated by the genetic algorithm in terms of quality of the
sampling (avoid steric clashes, good vina score, good smoothness).
"""

# Python
from __future__ import print_function
import random
import logging
import numpy as np
from math import sqrt
import copy
from zipfile import ZipFile, ZIP_STORED
import os
import pprint

# Chimera
import chimera
import Matrix as M
from Measure import inertia
from Molecule import atom_positions
from chimera import Xform as X
from FitMap.search import random_rotation

# GPATH
from gpath.genes import GeneProvider
from gpath import parse
from gpath.exceptions import LigandDirectoryAndTorsion, ProteinDirectoryAndNM

IDENTITY = ((1.0, 0.0, 0.0, 0.0),
            (0.0, 1.0, 0.0, 0.0),
            (0.0, 0.0, 1.0, 0.0))

ZERO = chimera.Point(0.0, 0.0, 0.0)

pp = pprint.PrettyPrinter(4)
logger = logging.getLogger(__name__)

def enable(**kwargs):
    kwargs = Pathway.validate(kwargs)
    return Pathway(**kwargs)

class Pathway(GeneProvider):
    """
    Pathway class

    Parameters
    ----------
    ligand : str
        Name of the GPathFinder ``molecule`` gene containing the ligand
        molecule.
    protein : str
        Name of the GPathFinder ``molecule`` gene containing the receptor
        molecule.
    torsion_gene : str, optional
        Name of the GPathFinder ``path_torsion`` gene to allow flexibility of
        the ligand.
    rotamers_gene : str, optional
        Name of the GPathFinder ``path_rotamers`` gene to allow conformational
        changes on the protein side-chain.
    radius_rotamers: float, optional, defaults to 3.0
        Maximum distance from any point of the ligand in every frame
        that is searched for possible rotamers  of the protein side-chain.
    nm_gene : str, optional
        Name of the GPathFinder ``path_normalmodes`` gene to allow the
        exploration of the protein global folding through normal modes analysis.
    origin : 3-item list or tuple of float, optional
        Coordinates to a specific geometric center at which the origin
        of the ligand will be set, if they are different from the actual
        position of the ligand in the .mol2 file.
    inertia_axes : max of 6-item list or tuple of int, optional
        Starting point of the ligand when generating binding pathways
        will be set at one of the indicated ends of the inertia axes of
        the protein.
    destination : 3-item list or tuple of float, optional
        It indicates the coordinates to the geometric binding site or the
        expected end point of the pathway.
        If not set by the user, GPathFinder will assume an unbinding scenario,
        calculating automatically the distance from the ligand to outside
        the protein.
    max_step_separation : float, optional
        Maximum distance in Angstroms from one point of the ligand in
        the pathway to the next one. If not set by the user, GPathFinder
        calculates the value from the size of the ligand.
    min_step_increment : float, optional
        Minimum distance increment in Angstroms from the ligand's origin
        that has to be the ligand in one frame of the pathway with
        respect of the ligand's distance from the origin of the previous
        frame of the pathway. If not set by the user, GPathFinder calculates
        the value as 1/5 by max_step_separation.
    frontier_margin : float, optional, defaults to 0.0
        Safety distance from the frontier of the protein to consider that
        the ligand is outside. Specially useful when using large ligands
        that can be stucked at the frontier.
    mut_pos_pb : float, optional, defaults to 0.10
        When a mutation occurs, this value is the probability of such
        mutation to be of the type `positions`, that is, the mutation
        changes the actual trajectory of the pathway. Warning: parameter
        for advanced users, usually the default value is correct for the
        vast majority of the systems.

    Attributes
    ----------
    allele : dict
        allele['positions'] = list of translation matrices of the ligand
                              along the frames of the pathway, as
                              explained in Notes.
        allele['rotations'] = list of rotation materices of the ligand
                              along the frames of the pathway, as
                              explained in Notes.
        allele['torsions'] = list of torsion alleles applied to the
                             ligand along the frames of the pathway (see
                             allele of GPathFinder ``path_torsion`` gene).
        allele['torsion_anchor'] = [str, int] ligand gene name and
                                   serial number of the anchor atom used
                                   in the torsion gene.
        allele['rotamers'] = list of data of the rotamers applied on the
                             protein along the frames of the pathway,
                             as explained in Notes.
        allele['normal_modes'] = list of normal modes alleles applied to
                                 the protein along the frames of the
                                 pathway (see allele of GPathFinder
                                 ``path_normalmodes`` gene).
    scores : list of dict
        Scores reported for every frame of the pathway. Each key of the
        dict correspond to one score (e.g. clashes).

    Notes
    -----
    ** Priority and configuration of the origin point of the ligand **

    The origin point/s of the ligand can be provided by three different
    manners:

        1. The ends of the inertia axes of the protein provided with the
           ``inertia_axes``.
        2. The coordinates set in ``origin`` parameter.
        3. The actual ubication of the ligand in its .mol2 file.

    The priority in case of conflict will be the order listed before.
    That is, if ``ìnertia_axes`` parameter is provided, this will
    override a possible ``origin`` parameter and the actual ubication in
    the .mol2 file.

    ** Details of the allele structure **

    - ``positions`` of the ligand along the pathway will be defined by
      translation matrices with the following shape:

        (
        (1, 0, 0, Tx),
        (0, 1, 0, Ty),
        (0, 0, 1, Tz)
        )

    - ``rotations`` of the ligand along the pathway will be defined by
      rotation matrices with the following shape:

        (
        (R1, R2, R3, 0),
        (R4, R5, R6, 0),
        (R7, R8, R9, 0)
        )

    - ``rotamers`` applied to the ligand at one frame of the pathway
      will be defined as 3 lists, each one containing:

        1. Residue numbers over which the rotamers are applied.
        2. Actual rotamer allele (see allele of GPathFinder ``path_rotamers`` gene)
        3. Chi angles corresponding to each rotamer.
    """

    _validate = {
        parse.Required('ligand'): parse.Molecule_name,
        parse.Required('protein'): parse.Molecule_name,
        'torsion_gene': parse.Molecule_name,
        'rotamers_gene': parse.Molecule_name,
        'radius_rotamers': parse.All(parse.Coerce(float), parse.Range(min=0)),
        'nm_gene': parse.Molecule_name,
        'origin': parse.Coordinates,
        'inertia_axes': [parse.All(parse.Coerce(int),
                            parse.Range(min=0, max=5))],
        'destination': parse.Coordinates,
        'max_step_separation': parse.Coerce(float),
        'min_step_increment': parse.Coerce(float),
        'frontier_margin': parse.Coerce(float),
        'mut_pos_pb': parse.Coerce(float),
        }

    def __init__(self, ligand, protein, torsion_gene=None, rotamers_gene=None,
                radius_rotamers=3.0, nm_gene=None,  origin=None,
                inertia_axes=None, destination=None, max_step_separation=None,
                min_step_increment=None, frontier_margin=0.0, mut_pos_pb=0.1, **kwargs):
        GeneProvider.__init__(self, **kwargs)
        self.ligand = ligand
        self.protein = protein
        self.torsion_gene = torsion_gene
        self.rotamers_gene = rotamers_gene
        self.radius_rotamers = radius_rotamers
        self.nm_gene = nm_gene
        self.destination = destination
        self.origin = origin if origin else None
        self.inertia_axes = inertia_axes if inertia_axes else None
        self.max_separation = max_step_separation
        self.min_increment = min_step_increment
        self.frontier_margin = frontier_margin
        self.mut_pos_pb = mut_pos_pb
        self.act_rotamers = []
        self.which_evaluations = []

    def __ready__(self):
        # Check if ligand conformers and torsion are used simultaneously
        if len(self.ligand_g.catalog) > 1 and self.torsion_gene:
            raise LigandDirectoryAndTorsion('If you set the Ligand molecule to be a directory, you can not use a path_torsion gene simultaneously')
        # Check if protein conformers and NM are used simultaneously
        if len(self.protein_g.catalog) > 1 and self.nm_gene:
            raise ProteinDirectoryAndNM('If you set the Protein molecule to be a directory, you can not use a path_normalmodes gene simultaneously')

        # Calculate inertia axes of the protein for further use
        self.p_axes, self.p_d2, self.p_center = inertia.atoms_inertia(self.protein_mol.atoms)
        self.p_elen = [a + self.frontier_margin for a in inertia.inertia_ellipsoid_size(self.p_d2)]
        self.p_axes = np.array(self.p_axes)

        if self.inertia_axes: #Prepare origin from inertia axes
            origin_points = []
            for axis, length in zip(self.p_axes, self.p_elen):
                origin_points.append((axis*length)+self.p_center)
                origin_points.append((-axis*length)+self.p_center)
            points = []
            for i, point in enumerate(origin_points):
                if i in self.inertia_axes:
                    points.append(point)
            self.origin = random.choice(points)
        elif not self.origin:
            #Origin is the geometric center of the ligand if not
            #provided explicitly
            self.origin = self.ligand_center

        #Calculate min_increment and max_separation if not set by user
        if not self.max_separation:
            axes, d2, center = inertia.atoms_inertia(self.ligand_mol.atoms)
            elen = [a for a in inertia.inertia_ellipsoid_size(d2)]
            self.max_separation = min(elen) if min(elen) > 1.0 else 1.0
        if self.min_increment is None:
            self.min_increment = 1.0 * self.max_separation / 5.0

        #Operations permitted to improve the path in mutate
        self._op_roulette = ['rotation']
        if self.torsion_gene:
            self._op_roulette.append('torsion')
        if self.rotamers_gene:
            self._op_roulette.append('rotamers')
        if self.nm_gene:
            self._op_roulette.append('nm')
        if len(self.protein_g.catalog) > 1:
            self._op_roulette.append('protein') #Conformers for the protein
        if len(self.ligand_g.catalog) > 1:
            self._op_roulette.append('ligand') #Conformers for the ligand

        #Making the initial allele
        self.allele = {}
        self.allele['positions'] = [((1.0, 0.0, 0.0, self.origin[0]),
                                    (0.0, 1.0, 0.0, self.origin[1]),
                                    (0.0, 0.0, 1.0, self.origin[2]))]
        self.allele['rotations'] = [IDENTITY]
        self.allele['min_increment'] = self.min_increment
        self.allele['max_separation'] = self.max_separation
        if self.torsion_gene:
            self.allele['torsions'] = [[]]
            self.allele['torsion_anchor'] = list(self.torsion_g._anchor)
        if self.rotamers_gene:
            self.allele['rotamers'] = [ [[],[],[]] ]
        if self.nm_gene:
            self.allele['normal_modes'] = [None]
        self.allele['protein'] = [0] #Conformer for the protein
        self.allele['ligand'] = [0] #Conformer for the ligand
        self.allele['coord_residues'] = [[]] #To store possible coord residues when evaluating metal_sites
        self.allele['mate_torsions'] = 0  #To register if mate is being useful
        self.allele['mate_positions'] = 0 #To register if mate is being useful
        self.allele['mate_rotations'] = 0 #To register if mate is being useful
        self.allele['mutate_positions'] = 0 #To register if mutate is being useful
        self.allele['mutate_torsions'] = 0 #To register if mutate is being useful
        self.allele['mutate_rotations'] = 0 #To register if mutate is being useful
        self.allele['mutate_rotamers'] = 0 #To register if mutate is being useful
        self.allele['mutate_nm'] = 0 #To register if mutate is being useful
        self.allele['mutate_protein'] = 0 #To register if mutate is being useful
        self.allele['mutate_ligand'] = 0 #To register if mutate is being useful
        self.scores = [{}] #To store scores at the evaluation stage

        #Initial pathway
        last_point = [x[3] for x in self.allele['positions'][-1]]
        if self.destination:
            not_arrived = distance(self.destination, last_point) > (self.max_separation)
        else:
            y = np.array([a-b for a,b in zip(last_point, self.p_center)])
            y = np.absolute(self.p_axes.dot(y))
            not_arrived = is_inside_ellipsoid(self.p_elen, y)

        while not_arrived:
            new_point = random_extension_point(last_point, self.max_separation,
                                self.origin, min_increment=self.min_increment,
                                destination=self.destination)
            if new_point is not None:
                self.allele['positions'].append(new_point)
                self.allele['rotations'].append(random_rotation())
                if self.torsion_gene:
                    self.torsion_g.gp_mutate(1.0)
                    self.allele['torsions'].append(copy.deepcopy(self.torsion_g.allele))
                if self.nm_gene:
                    self.allele['normal_modes'].append(random.randint(0,
                                    self.parent.genes[self.nm_gene].n_samples))
                self.allele['protein'].append(random.randint(0, len(self.protein_g.catalog)-1))
                self.allele['ligand'].append(random.randint(0, len(self.ligand_g.catalog)-1))
                if self.rotamers_gene:
                    self.allele['rotamers'].append([[],[],[]])
                self.allele['coord_residues'].append([])
                self.scores.append({})
                last_point = [x[3] for x in self.allele['positions'][-1]]
            if self.destination:
                not_arrived = distance(self.destination, last_point) > (self.max_separation)
            else:
                y = np.array([a-b for a,b in zip(last_point, self.p_center)])
                y = np.absolute(self.p_axes.dot(y))
                not_arrived = is_inside_ellipsoid(self.p_elen, y)

    def express(self):
        """
        Expression is controlled in gp_express to avoid unnecessary
        calls and supervise what frames are expressed.
        """
        pass

    def unexpress(self):
        """
        Unexpression is controlled in gp_unexpress to avoid unnecessary
        calls and supervise what frames are expressed.
        """
        pass

    def gp_express(self, i, with_rotamers=True, smoothness=False):
        """
        For the demanded frame ``i``:
            1. Apply translation, rotation and torsions to the ligand.
            2. Apply, if present, rotamers and normal modes to the
               protein.
        Optional parameter ``ẁith_rotamers`` controls if rotamers have
        to be expressed or not (useful when a rotamer actualization is
        pending).
        Optional parameter ``smoothness`` is used to express only the
        ligand transformations in order to calculate RMSD between frames.
        If True, only torsions and rotations of the ligand are applied.
        """
        ligand_sample = self.allele['ligand'][i]
        self.ligand_g.allele = self.ligand_g.catalog[ligand_sample]
        self.ligand_g._need_express = True
        self.ligand_g.express()

        protein_sample = self.allele['protein'][i]
        self.protein_g.allele = self.protein_g.catalog[protein_sample]
        self.protein_g._need_express = True
        self.protein_g.express()

        or_x, or_y, or_z = self.ligand_center
        to_zero = ((1.0, 0.0, 0.0, -or_x),
                    (0.0, 1.0, 0.0, -or_y),
                    (0.0, 0.0, 1.0, -or_z))

        rotations = self.allele['rotations']
        if smoothness:
            matrices = (IDENTITY,) + (rotations[i],) + (to_zero,)
        else:
            positions = self.allele['positions']
            matrices = (positions[i],) + (rotations[i],) + (to_zero,)
        matrices = M.multiply_matrices(*matrices)
        self.ligand_mol.openState.xform = M.chimera_xform(matrices)

        if self.torsion_gene:
            self.torsion_g.allele = copy.deepcopy(self.allele['torsions'][i])
            self.torsion_g._need_express = True
            self.torsion_g.express()

        if self.nm_gene and not smoothness:
            sample_number = self.allele['normal_modes'][i]
            if sample_number:
                self.nm_g.allele = self.nm_g.NORMAL_MODES_SAMPLES[sample_number]
                self.nm_g._need_express = True
                self.nm_g.express()

        if self.rotamers_gene and with_rotamers and not smoothness:
            if i in self.act_rotamers:
                #Actualize rotamers
                self.actualize_rotamers(i)

            if self.allele['rotamers'][i][1]:
                residues = []
                for rot in self.allele['rotamers'][i][0]:
                    a, b = rot.split('/')
                    residues.append([a, int(b)])

                self.rotamers_g._residues = residues
                self.rotamers_g.__ready__()
                self.rotamers_g.allele = self.allele['rotamers'][i][1]
                for (molname, pos), residue in self.rotamers_g.residues.items():
                    self.rotamers_g.allele[residues.index([molname, pos])] = self.allele['rotamers'][i][1][residues.index([molname, pos])]
                self.rotamers_g._need_express = True
                self.rotamers_g.express()

    def gp_unexpress(self, i, smoothness=False):
        """
        For the demanded frame ``i``:
            1. Undo translation, rotation and torsions made to the
               ligand.
            2. Undo, if present, rotamers and normal modes changes made
               to the protein.
        Optional parameter ``smoothness`` is used to unexpress only the
        ligand trasnformations when used to calculate RMSD of the ligand
        between frames. If True, only torsions and rotations of the
        ligand are unapplied.
        """
        self.ligand_mol.openState.xform = X()

        if self.torsion_gene:
            self.torsion_g.unexpress()
            self.torsion_g._need_express = False

        if self.rotamers_gene and not smoothness:
            if self.allele['rotamers'][i][0]:
                self.rotamers_g.unexpress()
                self.rotamers_g._need_express = False

        if self.nm_gene and not smoothness:
            sample_number = self.allele['normal_modes'][i]
            if sample_number:
                self.nm_g.unexpress()
                self.nm_g._need_express = False

        self.protein_g.unexpress()
        self.protein_g._need_express = False

        self.ligand_g.unexpress()
        self.ligand_g._need_express = False

    def mate(self, mate):
        """
        Randomnly choose an operation between [translations, rotations,
        torsions] of the ligand and make a mate of the allele of the two
        individuals.
        """
        #frame for self is the selected by score probability
        k = random.choice(self.scores[0].keys())
        scores_list = [sc[k] for sc in self.scores]
        m = min(scores_list)
        if m < 0:
            scores_list = [sc - m for sc in scores_list]
        scores_list = [sc + 1.0 for sc in scores_list] #to avoid divisions by 0 when calculating probs
        probs = [float(sc)/sum(scores_list[1:]) for sc in scores_list[1:]]
        i = np.random.choice(range(1, len(self.allele['positions'])), p=probs)

        #frame for mate is the nearest to the selected for self
        self_pos = [x[3] for x in self.allele['positions'][i]]
        dists = []
        for frame in range(1, len(mate.allele['positions'])):
            mate_pos = [x[3] for x in mate.allele['positions'][frame]]
            dists.append(distance(self_pos, mate_pos))
        j = dists.index(min(dists)) + 1 #dists doesn't contain distance for the frame 0

        if self.torsion_gene:
            op = random.choice([self.mate_torsions, self.mate_rotations])
        else:
            op = self.mate_rotations

        op(mate, i, j)

    def mate_rotations(self, mate, i, j):
        """
        For the frame ``i`` of self individual and ``j`` of the mate
        individual, make a interpolation of the rotations of the two
        ligands and assign the resulting matrix to both.
        """
        self.allele['mate_rotations'] = self.allele['mate_rotations'] + 1  #For register if mate is being useful
        mate.allele['mate_rotations'] = mate.allele['mate_rotations'] + 1  #For register if mate is being useful

        interp = interpolate_rotations(self.allele['rotations'][i],
                                            mate.allele['rotations'][j])
        self.allele['rotations'][i] = copy.deepcopy(interp)
        mate.allele['rotations'][j] = copy.deepcopy(interp)
        #Delete previous calculated scores
        self.scores[i] = {}
        mate.scores[j] = {}
        #Force rotamers actualization
        self.act_rotamers.append(i)
        mate.act_rotamers.append(j)

    def mate_torsions(self, mate, i, j):
        """
        For the frame ``i`` of self individual and ``j`` of the mate
        individual, make a torsion gp_mate (see ``path_torsion`` gene).
        """
        self.allele['mate_torsions'] = self.allele['mate_torsions'] + 1  #For register if mate is being useful
        mate.allele['mate_torsions'] = mate.allele['mate_torsions'] + 1  #For register if mate is being useful

        self.torsion_g.allele = copy.deepcopy(self.allele['torsions'][i])
        mate.torsion_g.allele = copy.deepcopy(mate.allele['torsions'][j])
        self.torsion_g.gp_mate(mate.torsion_g)
        self.allele['torsions'][i] = copy.deepcopy(self.torsion_g.allele)
        mate.allele['torsions'][j] = copy.deepcopy(mate.torsion_g.allele)
        #Delete previous calculated scores
        self.scores[i] = {}
        mate.scores[j] = {}
        #Force rotamers actualization
        self.act_rotamers.append(i)
        mate.act_rotamers.append(j)

    def mutate(self, indpb):
        """
        There are two possible mutations:
            1. Cut the pathway and remake  all the frames from the
               selected cutting frame (with probability mut_pos_pb).
            2. Change the [rotation, torsions, rotamers, normal modes]
               of one frame.
        The frame affected by the mutation is selected randomnly, with
        the probabilities of each frame weighted by their scores (more
        probability to be mutated if the frame has a worse score).
        The operation between the available pool is chosen randomly .
        """
        if random.random() < self.indpb:
            k = random.choice(self.scores[0].keys())
            scores_list = [sc[k] for sc in self.scores]
            m = min(scores_list)
            if m < 0:
                scores_list = [sc - m for sc in scores_list]
            scores_list = [sc + 1.0 for sc in scores_list] #to avoid divisions by 0 when calculating probs
            probs = [float(sc)/sum(scores_list[1:]) for sc in scores_list[1:]]
            i = np.random.choice(range(1, len(self.allele['positions'])), p=probs)

            n = random.random()
            if n < self.mut_pos_pb: #Modify the positions of the pathway
                i = np.random.choice(range(1, len(self.allele['positions']))) # No prob biass in this case
                self.allele['mutate_positions'] = self.allele['mutate_positions'] + 1 #For register if mutate is being useful
                #Delete frames from i
                self.allele['positions'] = self.allele['positions'][:i]
                self.allele['rotations'] = self.allele['rotations'][:i]
                if self.torsion_gene:
                    self.allele['torsions'] = self.allele['torsions'][:i]
                if self.rotamers_gene:
                    self.allele['rotamers'] = self.allele['rotamers'][:i]
                if self.nm_gene:
                    self.allele['normal_modes'] = self.allele['normal_modes'][:i]
                self.allele['protein'] = self.allele['protein'][:i]
                self.allele['ligand'] = self.allele['ligand'][:i]
                self.allele['coord_residues'] = self.allele['coord_residues'][:i]
                self.scores = self.scores[:i]
                #Generate new frames
                last_point = [x[3] for x in self.allele['positions'][-1]]
                if self.destination:
                    not_arrived = distance(self.destination, last_point) > (self.max_separation)
                else:
                    y = np.array([a-b for a,b in zip(last_point, self.p_center)])
                    y = np.absolute(self.p_axes.dot(y))
                    not_arrived = is_inside_ellipsoid(self.p_elen, y)

                while not_arrived:
                    new_point = random_extension_point(last_point, self.max_separation,
                                        self.origin, min_increment=self.min_increment,
                                        destination=self.destination)
                    if new_point is not None:
                        self.allele['positions'].append(new_point)
                        self.allele['rotations'].append(random_rotation())
                        if self.torsion_gene:
                            self.torsion_g.gp_mutate(1.0)
                            self.allele['torsions'].append(copy.deepcopy(self.torsion_g.allele))
                        if self.nm_gene:
                            self.allele['normal_modes'].append(random.randint(0,
                                            self.parent.genes[self.nm_gene].n_samples))
                        self.allele['protein'].append(random.randint(0, len(self.protein_g.catalog)-1))
                        self.allele['ligand'].append(random.randint(0, len(self.ligand_g.catalog)-1))
                        if self.rotamers_gene:
                            self.allele['rotamers'].append([[],[],[]])
                        self.allele['coord_residues'].append([])
                        self.scores.append({})
                        last_point = [x[3] for x in self.allele['positions'][-1]]
                    if self.destination:
                        not_arrived = distance(self.destination, last_point) > (self.max_separation)
                    else:
                        y = np.array([a-b for a,b in zip(last_point, self.p_center)])
                        y = np.absolute(self.p_axes.dot(y))
                        not_arrived = is_inside_ellipsoid(self.p_elen, y)
            else:
                #Try to modify rotation/torsion/rotamers/nm/sample protein
                operation = random.choice(self._op_roulette)
                if operation == 'rotation':
                    self.allele['mutate_rotations'] = self.allele['mutate_rotations'] + 1 #For register if mutate is being useful
                    self.allele['rotations'][i] = random_rotation()
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'torsion':
                    self.allele['mutate_torsions'] = self.allele['mutate_torsions'] + 1 #For register if mutate is being useful
                    self.torsion_g.gp_mutate(1.0)
                    self.allele['torsions'][i] = copy.deepcopy(self.torsion_g.allele)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'rotamers':
                    self.allele['mutate_rotamers'] = self.allele['mutate_rotamers'] + 1 #For register if mutate is being useful
                    self.scores[i] = {} #put the score to 0 and allow path_scoring to recalculate rotamers
                     #Search which rotamers are inside the search radius
                    residues = set()
                    atoms = [a for a in surrounding_atoms(self.ligand_mol, self.protein_mol, self.radius_rotamers)]
                    for atom in atoms:
                        residues.add(atom.residue.id.position)
                    residues_list = ['{}/{}'.format(self.protein, r) for r in residues]
                    residues = []
                    for rot in residues_list:
                        a, b = rot.split('/')
                        residues.append([a, int(b)])
                    #Create the new rotamers allele and express it to store in path allele
                    self.rotamers_g._residues = residues
                    self.rotamers_g.__ready__()
                    self.rotamers_g._need_express = True
                    self.rotamers_g.express()
                    #Store the new rotamer data in pathway gene
                    self.allele['rotamers'][i][0] = []
                    self.allele['rotamers'][i][1] = []
                    self.allele['rotamers'][i][2] = []
                    for ((molname, pos), residue), a in \
                                                    zip(self.rotamers_g.residues.items(),
                                                    self.rotamers_g.allele):
                        res_name = str(molname)+'/'+str(pos)
                        self.allele['rotamers'][i][0].append(res_name)
                        self.allele['rotamers'][i][1].append(a)
                        self.allele['rotamers'][i][2].append(self.rotamers_g.all_chis(residue))
                elif operation == 'nm':
                    self.allele['mutate_nm'] = self.allele['mutate_nm'] + 1 #For register if mutate is being useful
                    self.allele['normal_modes'][i] = random.randint(0,
                                    self.parent.genes[self.nm_gene].n_samples)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'protein':
                    self.allele['mutate_protein'] = self.allele['mutate_protein'] + 1 #For register if mutate is being useful
                    self.allele['protein'][i] = random.randint(0, len(self.protein_g.catalog)-1)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)
                elif operation == 'ligand':
                    self.allele['mutate_ligand'] = self.allele['mutate_ligand'] + 1 #For register if mutate is being useful
                    self.allele['ligand'][i] = random.randint(0, len(self.ligand_g.catalog)-1)
                    self.scores[i] = {}
                    self.act_rotamers.append(i)

    def actualize_rotamers(self, i):
        """
        Actualize the rotamers list and allele when a position/
        orientation/torsions of a ligand have changed (it could be
        possible that with the new ligand expression the surrounding
        rotamers would be different, so it is necessary to actualize the
        information). Old rotamer allele is maintained when possible
        (i.e. the rotamers that are present in the previous and the new
        sourrondings of the ligand).
        """
        #Search which rotamers are inside the search radius
        residues = set()
        atoms = [a for a in surrounding_atoms(self.ligand_mol, self.protein_mol, self.radius_rotamers)]
        for atom in atoms:
            residues.add(atom.residue.id.position)
        residues_list = ['{}/{}'.format(self.protein, r) for r in residues]
        residues = []
        for rot in residues_list:
            a, b = rot.split('/')
            residues.append([a, int(b)])

        #Create the new rotamers allele and express i
        self.rotamers_g._residues = residues
        self.rotamers_g.__ready__()
        self.rotamers_g._need_express = True
        self.rotamers_g.express()

        #Store the new rotamer data in pathway gene (only new residues)
        old_res = [[],[],[]]
        old_res[0] = list(self.allele['rotamers'][i][0])
        old_res[1] = list(self.allele['rotamers'][i][1])
        old_res[2] = list(self.allele['rotamers'][i][2])
        self.allele['rotamers'][i][0] = []
        self.allele['rotamers'][i][1] = []
        self.allele['rotamers'][i][2] = []
        new_res = [str(molname)+'/'+str(pos) for ((molname, pos), residue), a in \
                    zip(self.rotamers_g.residues.items(),
                        self.rotamers_g.allele)]
        #Incorporate old residues present in the new region as is
        for j, res_name in enumerate(old_res[0]):
            if res_name in new_res:
                self.allele['rotamers'][i][0].append(res_name)
                self.allele['rotamers'][i][1].append(old_res[1][j])
                self.allele['rotamers'][i][2].append(old_res[2][j])
        #Incorporate new residues of the new region
        for ((molname, pos), residue), a in \
                    zip(self.rotamers_g.residues.items(),
                        self.rotamers_g.allele):
            res_name = str(molname)+'/'+str(pos)
            if res_name not in self.allele['rotamers'][i][0]:
                self.allele['rotamers'][i][0].append(res_name)
                self.allele['rotamers'][i][1].append(a)
                self.allele['rotamers'][i][2].append(self.rotamers_g.all_chis(residue))

        self.rotamers_g.unexpress()
        self.rotamers_g._need_express = False

        self.act_rotamers.remove(i) #Rotamers are actualized

    def write(self, path, name, *args, **kwargs):
        """
        Each individual output is a .zip file that contains:
            1. A .txt file with the allele data (allele.txt).
            2. A .txt file with the scoring data (scores.txt).
            3. One .pdb file for each frame of the calculated pathway
               (named frame_xxx.pdb, where xxx is the number of frame).
        """
        fullname = os.path.join(path, '{}_{}.zip'.format(name, self.name))
        fullname = self.create_output_zip(path, fullname)
        return fullname

    def create_output_zip(self, path, fullname):
        #Txt file with the allele
        allelename = os.path.join(path, 'allele.txt')
        with open(allelename, 'w') as f:
            f.write(pp.pformat(self.allele))

        #Txt file with the scores
        scoresname = os.path.join(path, 'scores.txt')
        with open(scoresname, 'w') as f:
            f.write(pp.pformat(self.scores))

        #Pdb file with the trajectory of the pathway
        trajname = os.path.join(path, 'trajectory.pdb')
        points = self.allele['positions']
        mol = chimera.Molecule()
        r = mol.newResidue("path", " ", 1, " ")
        atoms = []
        atom = mol.newAtom("0", chimera.Element("H"))
        atoms.append(atom)
        point = [x[3] for x in points[0]]
        atoms[-1].setCoord(chimera.Point(*point))
        r.addAtom(atoms[-1])
        for i in range(1, len(points)):
            point = [x[3] for x in points[i]]
            atom = mol.newAtom(str(i), chimera.Element("H"))
            atoms.append(atom)
            atoms[-1].setCoord(chimera.Point(*point))
            r.addAtom(atoms[-1])
            mol.newBond(atoms[i-1], atoms[i])
        chimera.pdbWrite([mol], chimera.Xform(), trajname)

        with ZipFile(fullname, 'w', ZIP_STORED) as z:
            z.write(allelename, os.path.basename(allelename))
            os.remove(allelename)
            z.write(scoresname, os.path.basename(scoresname))
            os.remove(scoresname)
            z.write(trajname, os.path.basename(trajname))
            os.remove(trajname)

            #Pdb files with the actual frames of the pathway
            for i in range(0, len(self.allele['positions'])):

                self.gp_express(i)
                framename = os.path.join(path, "frame_{:03d}.pdb".format(i))
                chimera.pdbWrite([self.protein_mol, self.ligand_mol],
                                    chimera.Xform(), framename)

                # Prepare pdb files with format required by Chimera MD movie import
                with open(framename) as frame_file:
                    pdb_file_lines = []
                    line = frame_file.readline()
                    while line:
                        if line[:4] == "ATOM" or line[:6] == "HETATM":
                            pdb_file_lines.append(line)
                        line = frame_file.readline()

                os.remove(framename)
                with open(framename, 'w') as f:
                    prev_atom = 0
                    prev_res = 0
                    last_change_res = 0
                    for line in pdb_file_lines:
                        current_atom, current_res = int(line[6:11]), int(line[22:26])
                        if current_atom > prev_atom:
                            prev_atom = current_atom
                        else:
                            current_atom = prev_atom + 1
                            prev_atom = current_atom
                        if current_res > prev_res:
                            prev_res = current_res
                            last_change_res = current_res
                        elif current_res != last_change_res:
                            last_change_res = current_res
                            prev_res += 1
                            current_res = prev_res
                        else:
                            current_res = prev_res
                        line = line[:6] + '{:5d}'.format(current_atom) + line[11:22] + '{:4d}'.format(current_res) + line[26:]
                        f.write("%s" % line)

                z.write(framename, os.path.basename(framename))
                os.remove(framename)
                self.gp_unexpress(i)

        return fullname

    #####
    @property
    def ligand_mol(self):
        return self.parent.find_molecule(self.ligand).compound.mol

    @property
    def protein_mol(self):
        return self.parent.find_molecule(self.protein).compound.mol

    @property
    def ligand_g(self):
        return self.parent.find_molecule(self.ligand)

    @property
    def protein_g(self):
        return self.parent.find_molecule(self.protein)

    @property
    def torsion_g(self):
        if self.torsion_gene:
            return self.parent.genes[self.torsion_gene]
        else:
            return None

    @property
    def rotamers_g(self):
        if self.rotamers_gene:
            return self.parent.genes[self.rotamers_gene]
        else:
            return None

    @property
    def nm_g(self):
        if self.nm_gene:
            return self.parent.genes[self.nm_gene]
        else:
            return None

    @property
    #Geometric center of the ligand
    def ligand_center(self):
        coordinates = atom_positions(self.ligand_mol.atoms)
        c = np.average(coordinates, axis=0)
        return c


#############
# Some useful functions
def random_extension_point(center, r, origin, min_increment=0, destination=None):
    """
    Calculate a new random coordinates that acomplishes the following
    conditions:
        1. The distance from the ``center`` has to be <= r.
        2. The distance from the origin has to be >= than that from the
           ``center`` to origin plus ``min_increment`` when calculating
           pathways without a ``destination``.
        3. If a ``destination`` is set, the new point has to be nearer
           to this destination than the ``center`` point.

    Parameters
    ----------
    center : 3-tuple of float
        Coordinates of the previous point (origin of the search sphere).
    r : float
        Radius of the search sphere (maximum distance at will the new
        point will be from ``center``).
    origin : 3-tuple of float
        Coordinates of the origin point of the pathway.
    min_increment : float, optional, defaults to 0
        Minimum increment of distance from the origin / to the destination
        of the new point respect to th previous one.
    destination : float or 3-tuple of float, optional
        Distance from the origin or coordinates of the destination point
        of the pathway.

    Returns
    -------
    A 3-tuple of float representing the coordinates that accomplish all
    the criteria. If a valid point is not obtained in 100 iterations,
    returns None.

    Notes
    -----
    The combination of the ``r`` (max_step_separation) and
    ``min_increment`` will control the distance of the different frames
    of the pathway between them and the possibility of deviation from the
    straight line (as higher the difference between these parameters,
    more deviation from the straight line will be allowed).
    """
    for i in range(0, 100):
        new_point = [random.uniform(-r, r) for m in center]
        new_point = [x + y for x, y in zip(new_point, center)]
        x, y, z = new_point
        if not isinstance(destination, list):
            if ((distance(new_point, origin) >= (distance(center, origin) +
                                                min_increment))
                and (distance(new_point, center) <= r)):
                return ((1.0, 0.0, 0.0, x),
                        (0.0, 1.0, 0.0, y),
                        (0.0, 0.0, 1.0, z))
        else:
            if (((distance(new_point, destination) +
                    min_increment) < distance(center, destination))
                and (distance(new_point, center) <= r)):
                return ((1.0, 0.0, 0.0, x),
                        (0.0, 1.0, 0.0, y),
                        (0.0, 0.0, 1.0, z))
    return None

def distance(a, b):
    ax, ay, az = a
    bx, by, bz = b
    return(sqrt((bx-ax)*(bx-ax)+(by-ay)*(by-ay)+(bz-az)*(bz-az)))

def interpolate_rotations(a, b):
    xf1 = M.chimera_xform(a)
    xf2 = M.chimera_xform(b)
    interp = M.xform_matrix(M.interpolate_xforms(xf1, ZERO, xf2, 0.5))
    return interp

def is_inside_ellipsoid(elen, point):
    discr = (point[0]/elen[0])**2 + (point[1]/elen[1])**2 + (point[2]/elen[2])**2
    if discr < 1.0:
        return True
    else:
        return False

def surrounding_atoms(ligand, protein, radius_rotamers):
    """
    Get atoms in the search zone, based on the molecule and the
    radius.
    """
    z = chimera.selection.ItemizedSelection()
    z.add([a for a in ligand.atoms])
    z.merge(chimera.selection.REPLACE,
                    chimera.specifier.zone(z, 'atom', None,
                                            radius_rotamers,
                                            [ligand, protein]))
    atom_list = [at for at in z.atoms() if at not in ligand.atoms]
    return atom_list
#!/usr/bin/env python # -*- coding: utf-8 -*- ############## # GPathFinder: Identification of ligand pathways by a multi-objective # genetic algorithm # # https://github.com/insilichem/gpathfinder # # Copyright 2019 <NAME>, <NAME>, # <NAME>, <NAME>, # <NAME> and <NAME> # # 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. # See the License for the specific language governing permissions and # limitations under the License. ############## GPathFinder allows to generate approximate ligand binding/unbinding pathways to/from the binding site of a receptor (typically a protein).  This module must be used in conjunction with at least the following three genes: - A GPathFinder standard ``molecule`` gene for the ligand molecule. - A GPathFinder standard ``molecule`` gene for the protein molecule. - A GPathFinder ``path_torsion`` gene applied on the ligand to allow its flexibility.  Optionally, the following genes can be used as well: - A GPathFinder ``path_rotamers`` gene on the protein to allow conformational changes in its side-chain. - A GPathFinder``path_normal modes`` gene on the protein to allow the exploration of its global folding through normal modes analysis.  To assess the quality of the pathways generated by this pool of genes, the following objective must be used in the evaluation stage: - At least one ``path scoring`` objective to evaluate how good are the pathways generated by the genetic algorithm in terms of quality of the sampling (avoid steric clashes, good vina score, good smoothness). # Python # Chimera # GPATH Pathway class      Parameters     ----------     ligand : str         Name of the GPathFinder ``molecule`` gene containing the ligand         molecule.     protein : str         Name of the GPathFinder ``molecule`` gene containing the receptor         molecule.     torsion_gene : str, optional         Name of the GPathFinder ``path_torsion`` gene to allow flexibility of         the ligand.     rotamers_gene : str, optional         Name of the GPathFinder ``path_rotamers`` gene to allow conformational         changes on the protein side-chain.     radius_rotamers: float, optional, defaults to 3.0         Maximum distance from any point of the ligand in every frame         that is searched for possible rotamers  of the protein side-chain.     nm_gene : str, optional         Name of the GPathFinder ``path_normalmodes`` gene to allow the         exploration of the protein global folding through normal modes analysis.     origin : 3-item list or tuple of float, optional         Coordinates to a specific geometric center at which the origin         of the ligand will be set, if they are different from the actual         position of the ligand in the .mol2 file.     inertia_axes : max of 6-item list or tuple of int, optional         Starting point of the ligand when generating binding pathways         will be set at one of the indicated ends of the inertia axes of         the protein.     destination : 3-item list or tuple of float, optional         It indicates the coordinates to the geometric binding site or the         expected end point of the pathway.         If not set by the user, GPathFinder will assume an unbinding scenario,         calculating automatically the distance from the ligand to outside         the protein.     max_step_separation : float, optional         Maximum distance in Angstroms from one point of the ligand in         the pathway to the next one. If not set by the user, GPathFinder         calculates the value from the size of the ligand.     min_step_increment : float, optional         Minimum distance increment in Angstroms from the ligand's origin         that has to be the ligand in one frame of the pathway with         respect of the ligand's distance from the origin of the previous         frame of the pathway. If not set by the user, GPathFinder calculates         the value as 1/5 by max_step_separation.     frontier_margin : float, optional, defaults to 0.0         Safety distance from the frontier of the protein to consider that         the ligand is outside. Specially useful when using large ligands         that can be stucked at the frontier.     mut_pos_pb : float, optional, defaults to 0.10         When a mutation occurs, this value is the probability of such         mutation to be of the type `positions`, that is, the mutation         changes the actual trajectory of the pathway. Warning: parameter         for advanced users, usually the default value is correct for the         vast majority of the systems.      Attributes     ----------     allele : dict         allele['positions'] = list of translation matrices of the ligand                               along the frames of the pathway, as                               explained in Notes.         allele['rotations'] = list of rotation materices of the ligand                               along the frames of the pathway, as                               explained in Notes.         allele['torsions'] = list of torsion alleles applied to the                              ligand along the frames of the pathway (see                              allele of GPathFinder ``path_torsion`` gene).         allele['torsion_anchor'] = [str, int] ligand gene name and                                    serial number of the anchor atom used                                    in the torsion gene.         allele['rotamers'] = list of data of the rotamers applied on the                              protein along the frames of the pathway,                              as explained in Notes.         allele['normal_modes'] = list of normal modes alleles applied to                                  the protein along the frames of the                                  pathway (see allele of GPathFinder                                  ``path_normalmodes`` gene).     scores : list of dict         Scores reported for every frame of the pathway. Each key of the         dict correspond to one score (e.g. clashes).      Notes     -----     ** Priority and configuration of the origin point of the ligand **      The origin point/s of the ligand can be provided by three different     manners:          1. The ends of the inertia axes of the protein provided with the            ``inertia_axes``.         2. The coordinates set in ``origin`` parameter.         3. The actual ubication of the ligand in its .mol2 file.      The priority in case of conflict will be the order listed before.     That is, if ``ìnertia_axes`` parameter is provided, this will     override a possible ``origin`` parameter and the actual ubication in     the .mol2 file.      ** Details of the allele structure **      - ``positions`` of the ligand along the pathway will be defined by       translation matrices with the following shape:          (         (1, 0, 0, Tx),         (0, 1, 0, Ty),         (0, 0, 1, Tz)         )      - ``rotations`` of the ligand along the pathway will be defined by       rotation matrices with the following shape:          (         (R1, R2, R3, 0),         (R4, R5, R6, 0),         (R7, R8, R9, 0)         )      - ``rotamers`` applied to the ligand at one frame of the pathway       will be defined as 3 lists, each one containing:          1. Residue numbers over which the rotamers are applied.         2. Actual rotamer allele (see allele of GPathFinder ``path_rotamers`` gene)         3. Chi angles corresponding to each rotamer. # Check if ligand conformers and torsion are used simultaneously # Check if protein conformers and NM are used simultaneously # Calculate inertia axes of the protein for further use #Prepare origin from inertia axes #Origin is the geometric center of the ligand if not #provided explicitly #Calculate min_increment and max_separation if not set by user #Operations permitted to improve the path in mutate #Conformers for the protein #Conformers for the ligand #Making the initial allele #Conformer for the protein #Conformer for the ligand #To store possible coord residues when evaluating metal_sites #To register if mate is being useful #To register if mate is being useful #To register if mate is being useful #To register if mutate is being useful #To register if mutate is being useful #To register if mutate is being useful #To register if mutate is being useful #To register if mutate is being useful #To register if mutate is being useful #To register if mutate is being useful #To store scores at the evaluation stage #Initial pathway Expression is controlled in gp_express to avoid unnecessary         calls and supervise what frames are expressed. Unexpression is controlled in gp_unexpress to avoid unnecessary         calls and supervise what frames are expressed. For the demanded frame ``i``:             1. Apply translation, rotation and torsions to the ligand.             2. Apply, if present, rotamers and normal modes to the                protein.         Optional parameter ``ẁith_rotamers`` controls if rotamers have         to be expressed or not (useful when a rotamer actualization is         pending).         Optional parameter ``smoothness`` is used to express only the         ligand transformations in order to calculate RMSD between frames.         If True, only torsions and rotations of the ligand are applied. #Actualize rotamers For the demanded frame ``i``:             1. Undo translation, rotation and torsions made to the                ligand.             2. Undo, if present, rotamers and normal modes changes made                to the protein.         Optional parameter ``smoothness`` is used to unexpress only the         ligand trasnformations when used to calculate RMSD of the ligand         between frames. If True, only torsions and rotations of the         ligand are unapplied. Randomnly choose an operation between [translations, rotations,         torsions] of the ligand and make a mate of the allele of the two         individuals. #frame for self is the selected by score probability #to avoid divisions by 0 when calculating probs #frame for mate is the nearest to the selected for self #dists doesn't contain distance for the frame 0 For the frame ``i`` of self individual and ``j`` of the mate         individual, make a interpolation of the rotations of the two         ligands and assign the resulting matrix to both. #For register if mate is being useful #For register if mate is being useful #Delete previous calculated scores #Force rotamers actualization For the frame ``i`` of self individual and ``j`` of the mate         individual, make a torsion gp_mate (see ``path_torsion`` gene). #For register if mate is being useful #For register if mate is being useful #Delete previous calculated scores #Force rotamers actualization There are two possible mutations:             1. Cut the pathway and remake  all the frames from the                selected cutting frame (with probability mut_pos_pb).             2. Change the [rotation, torsions, rotamers, normal modes]                of one frame.         The frame affected by the mutation is selected randomnly, with         the probabilities of each frame weighted by their scores (more         probability to be mutated if the frame has a worse score).         The operation between the available pool is chosen randomly . #to avoid divisions by 0 when calculating probs #Modify the positions of the pathway # No prob biass in this case #For register if mutate is being useful #Delete frames from i #Generate new frames #Try to modify rotation/torsion/rotamers/nm/sample protein #For register if mutate is being useful #For register if mutate is being useful #For register if mutate is being useful #put the score to 0 and allow path_scoring to recalculate rotamers #Search which rotamers are inside the search radius #Create the new rotamers allele and express it to store in path allele #Store the new rotamer data in pathway gene #For register if mutate is being useful #For register if mutate is being useful #For register if mutate is being useful Actualize the rotamers list and allele when a position/         orientation/torsions of a ligand have changed (it could be         possible that with the new ligand expression the surrounding         rotamers would be different, so it is necessary to actualize the         information). Old rotamer allele is maintained when possible         (i.e. the rotamers that are present in the previous and the new         sourrondings of the ligand). #Search which rotamers are inside the search radius #Create the new rotamers allele and express i #Store the new rotamer data in pathway gene (only new residues) #Incorporate old residues present in the new region as is #Incorporate new residues of the new region #Rotamers are actualized Each individual output is a .zip file that contains:             1. A .txt file with the allele data (allele.txt).             2. A .txt file with the scoring data (scores.txt).             3. One .pdb file for each frame of the calculated pathway                (named frame_xxx.pdb, where xxx is the number of frame). #Txt file with the allele #Txt file with the scores #Pdb file with the trajectory of the pathway #Pdb files with the actual frames of the pathway # Prepare pdb files with format required by Chimera MD movie import ##### #Geometric center of the ligand ############# # Some useful functions Calculate a new random coordinates that acomplishes the following     conditions:         1. The distance from the ``center`` has to be <= r.         2. The distance from the origin has to be >= than that from the            ``center`` to origin plus ``min_increment`` when calculating            pathways without a ``destination``.         3. If a ``destination`` is set, the new point has to be nearer            to this destination than the ``center`` point.      Parameters     ----------     center : 3-tuple of float         Coordinates of the previous point (origin of the search sphere).     r : float         Radius of the search sphere (maximum distance at will the new         point will be from ``center``).     origin : 3-tuple of float         Coordinates of the origin point of the pathway.     min_increment : float, optional, defaults to 0         Minimum increment of distance from the origin / to the destination         of the new point respect to th previous one.     destination : float or 3-tuple of float, optional         Distance from the origin or coordinates of the destination point         of the pathway.      Returns     -------     A 3-tuple of float representing the coordinates that accomplish all     the criteria. If a valid point is not obtained in 100 iterations,     returns None.      Notes     -----     The combination of the ``r`` (max_step_separation) and     ``min_increment`` will control the distance of the different frames     of the pathway between them and the possibility of deviation from the     straight line (as higher the difference between these parameters,     more deviation from the straight line will be allowed). Get atoms in the search zone, based on the molecule and the     radius.