works approximately the same as old, but with less functionality. tiny...

works approximately the same as old, but with less functionality. tiny differences between old and new, when running with same seed. They do however diverge more as the searches progress.

bfgslinesearch_zlim.py

+from __future__ import print_function
+
+# ******NOTICE***************
+# optimize.py module by Travis E. Oliphant
+#
+# You may copy and use this module as you see fit with no
+# guarantee implied provided you keep this notice in all copies.
+# *****END NOTICE************
+
+import time
+import numpy as np
+from numpy import eye, absolute, sqrt, isinf
+from ase.utils.linesearch import LineSearch
+from ase.optimize.bfgslinesearch import BFGSLineSearch
+from ase.constraints import FixAtoms
+from ase.utils import basestring
+from ase.data import covalent_radii
+from ase.ga.utilities import get_mic_distance
+
+# These have been copied from Numeric's MLab.py
+# I don't think they made the transition to scipy_core
+
+# Modified from scipy_optimize
+abs = absolute
+pymin = min
+pymax = max
+__version__ = '0.1'
+
+
+class BFGSLineSearch_zlim(BFGSLineSearch):
+    def __init__(self, atoms, pos_init=None, restart=None, logfile='-', maxstep=.2,
+                 trajectory=None, c1=0.23, c2=0.46, alpha=10.0, stpmax=50.0,
+                 master=None, force_consistent=None, zlim=None, dmax_cov=1.5, rk=None):
+        """
+        add maximum displacement of single atoms to BFGSLineSearch:
+
+        dmax_cov: maximum distance the atom is alowed to move from it's initial position.
+        in units of it's covalent distance.
+        """
+
+        self.rk = rk  # for testing
+        
+        self.zlim = zlim
+        
+        self.dmax_cov = dmax_cov
+        self.d_cov = 2*np.array([covalent_radii[num] for num in atoms.get_atomic_numbers()])
+        if pos_init is not None:
+            self.pos_init = pos_init
+        else:
+            self.pos_init = np.copy(atoms.positions)
+
+        self.cell = atoms.get_cell()
+        self.pbc = atoms.get_pbc()
+
+        BFGSLineSearch.__init__(self, atoms, restart=restart, logfile=logfile, maxstep=maxstep,
+                                trajectory=trajectory, c1=c1, c2=c2, alpha=alpha, stpmax=stpmax,
+                                master=master, force_consistent=force_consistent)
+
+    def converged(self, forces=None):
+        """Did the optimization converge?"""
+        if forces is None:
+            forces = self.atoms.get_forces()
+        if hasattr(self.atoms, 'get_curvature'):
+            return ((forces**2).sum(axis=1).max() < self.fmax**2 and
+                    self.atoms.get_curvature() < 0.0)
+
+        # Check if stop due to large displacement
+        if self.dmax_cov is not None:
+            #max_covDisplace = np.max(np.sqrt(((self.pos_init-self.atoms.positions)**2).sum(axis=1))/self.d_cov)
+            d = np.array([get_mic_distance(p1,p2,self.cell,self.pbc) for p1,p2 in zip(self.pos_init,self.atoms.get_positions())])
+            max_covDisplace = np.max(d/self.d_cov)
+            if max_covDisplace > self.dmax_cov:
+                print('RELAX STOPPED on rank {} - max relax dist (>{} covDist) reached.'.format(self.rk, self.dmax_cov), flush=True)
+                return True
+        # Check if stop due to zlim
+        if self.zlim is not None:
+            # get indices of non-fixed atoms
+            indices = np.arange(self.atoms.get_number_of_atoms())
+            for constraint in self.atoms.constraints:
+                if isinstance(constraint, FixAtoms):
+                    indices_fixed = constraint.get_indices()
+                    indices = np.delete(np.arange(self.atoms.get_number_of_atoms()), indices_fixed)
+            pos_z = self.atoms.positions[indices,2]
+            if np.any(pos_z < self.zlim[0]) or np.any(pos_z > self.zlim[1]):
+                print('RELAXATION STOPPED on rank {} due to - zlim crossed.'.format(self.rk))
+                return True
+        return (forces**2).sum(axis=1).max() < self.fmax**2
+

candidate_operations/candidateGenerator.py

+import numpy as np
+
+
+
+class candidateGenerator():
+    def __init__(self):
+        pass
+
+
+
+
+class candidateOperation():
+    def __init__(self):
+        pass

candidate_operations/offspring_creator.py

+"""Base module for all operators that create offspring."""
+import numpy as np
+from ase import Atoms
+
+
+class OffspringCreator(object):
+    """Base class for all procreation operators
+
+    Parameters:
+
+    verbose: Be verbose and print some stuff
+
+    """
+
+    def __init__(self, verbose=False, num_muts=1):
+        self.descriptor = 'OffspringCreator'
+        self.verbose = verbose
+        self.min_inputs = 0
+        self.num_muts = num_muts
+
+    def get_min_inputs(self):
+        """Returns the number of inputs required for a mutation,
+        this is to know how many candidates should be selected
+        from the population."""
+        return self.min_inputs
+
+    def get_new_individual(self, parents):
+        """Function that returns a new individual.
+        Overwrite in subclass."""
+        raise NotImplementedError
+
+    def finalize_individual(self, indi):
+        """Call this function just before returning the new individual"""
+        indi.info['key_value_pairs']['origin'] = self.descriptor
+
+        return indi
+
+    @classmethod
+    def initialize_individual(cls, parent, indi=None):
+        """Initializes a new individual that inherits some parameters
+        from the parent, and initializes the info dictionary.
+        If the new individual already has more structure it can be
+        supplied in the parameter indi."""
+        if indi is None:
+            indi = Atoms(pbc=parent.get_pbc(), cell=parent.get_cell())
+        else:
+            indi = indi.copy()
+        # key_value_pairs for numbers and strings
+        indi.info['key_value_pairs'] = {'extinct': 0}
+        # data for lists and the like
+        indi.info['data'] = {}
+
+        return indi
+
+
+class OperationSelector(object):
+    """Class used to randomly select a procreation operation
+    from a list of operations.
+
+    Parameters:
+
+    probabilities: A list of probabilities with which the different
+        mutations should be selected. The norm of this list
+        does not need to be 1.
+
+    oplist: The list of operations to select from.
+    """
+
+    def __init__(self, probabilities=None, oplist=None):
+        assert len(probabilities) == len(oplist)
+        self.oplist = oplist
+        self.rho = np.cumsum(probabilities)
+        
+    def __get_index__(self, extra=False):
+        v = np.random.random() * self.rho[-1]
+        for i in range(len(self.rho)):
+            if self.rho[i] > v:
+                return i
+
+    def get_new_individual(self, candidate_list, debugger=None):
+        """Choose operator and use it on the candidate. """
+        to_use = self.__get_index__()
+        a, mut = self.oplist[to_use].get_new_individual(candidate_list)
+        return a, mut
+
+    def get_operator(self):
+        """Choose operator and return it."""
+        to_use = self.__get_index__()
+        return self.oplist[to_use]

candidate_operations/standardmutations_mkb.py

+""" A collection of mutations that can be used. """
+
+import numpy as np
+from math import ceil, cos, sin, pi
+from ase.ga.utilities import atoms_too_close
+from ase.ga.utilities import atoms_too_close_two_sets
+from mutations.offspring_creator import OffspringCreator
+from ase import Atoms, Atom
+
+def pos_add_sphere(rattle_strength):
+    # Rattle within a sphere
+    r = rattle_strength * np.random.rand()**(1/3)
+    theta = np.random.uniform(low=0, high=2*np.pi)
+    phi = np.random.uniform(low=0, high=np.pi)
+    pos_add = r * np.array([np.cos(theta)*np.sin(phi),
+                            np.sin(theta)*np.sin(phi),
+                            np.cos(phi)])
+    return pos_add
+
+def pos_add_cylinder(rattle_strength, direction, atoms):
+    direction = int(direction)
+    cylinder_length = atoms.cell[direction][direction]
+    # Rattle within a cylinder
+    r = rattle_strength * np.random.rand()**(1/2)
+    theta = np.random.uniform(low=0, high=2*np.pi)
+
+    pos_add_axis = np.random.uniform(low=-0.5*cylinder_length, high=0.5*cylinder_length)
+    pos_add_radial = r * np.array([np.cos(theta),
+                                   np.sin(theta)])
+    if direction == 0:
+        pos_add = np.array([pos_add_axis,
+                            pos_add_radial[0],
+                            pos_add_radial[1]])
+    elif direction == 1:
+        pos_add = r * np.array([pos_add_radial[0],
+                                pos_add_axis,
+                                pos_add_radial[1]])
+    elif direction == 2:
+        pos_add = r * np.array([pos_add_radial[0],
+                                pos_add_radial[1],
+                                pos_add_axis])
+    else:
+        raise RuntimeError('cylinder direction can only be 0 or 1 or 2')
+    return pos_add
+
+
+class RattleMutation(OffspringCreator):
+    """ 
+    An implementation of the rattle mutation as described in
+    R.L. Johnston Dalton Transactions, Vol. 22,
+    No. 22. (2003), pp. 4193-4207
+
+    Parameters:
+
+    blmin: Dictionary defining the minimum distance between atoms
+    after the rattle.
+    n_top: Number of atoms optimized by the GA.
+    rattle_strength: Strength with which the atoms are moved.
+    rattle_prop: The probability with which each atom is rattled.
+    min_z: minimum allowed z-coordinate for rattled atoms.
+    rattle_mode: 'sphere' or 'cylinder'
+    """
+    def __init__(self, blmin, n_top, rattle_strength=0.8,
+                 rattle_prop=0.4, min_z=None, use_max_dist=True, verbose=False,
+                 rattle_mode='sphere', cluster=False, cylinder_direction=0,
+                 descriptor='RattleMutation'):
+        OffspringCreator.__init__(self, verbose)
+        self.blmin = blmin
+        self.n_top = n_top
+        self.rattle_strength = rattle_strength
+        self.rattle_prop = rattle_prop
+        self.min_z = min_z
+        self.use_max_dist = use_max_dist
+        self.rattle_mode = rattle_mode
+        self.cluster = cluster
+        self.cylinder_direction = cylinder_direction
+        self.descriptor = descriptor
+        self.min_inputs = 1
+
+    def get_new_individual(self, parents, **kwargs):
+        f = parents[0]
+
+        indi = self.mutate(f)
+        if indi is None:
+            return indi, 'mutation: rattle'
+            
+        indi = self.initialize_individual(f, indi)
+        #indi.info['data']['parents'] = [f.info['confid']]  # remove for GP-search
+
+        return self.finalize_individual(indi), 'mutation: rattle'
+
+    def mutate(self, atoms):
+        num = atoms.numbers
+        atoms_mutated = atoms.copy()
+        Natoms = len(atoms)
+        Nslab = Natoms - self.n_top
+        for i in range(Nslab, Natoms):
+            if np.random.random() < self.rattle_prop:
+                posi_0 = np.copy(atoms_mutated.positions[i])
+                tc = False
+                tf = False
+                for k in range(100):
+                    atoms_mutated.positions[i] = posi_0
+
+                    # Make random displacement of atom i
+                    if self.rattle_mode == 'sphere':
+                        pos_add = pos_add_sphere(self.rattle_strength)
+                    elif self.rattle_mode == 'cylinder':
+                        pos_add = pos_add_cylinder(self.rattle_strength, self.cylinder_direction, atoms)
+                    atoms_mutated.positions[i] += pos_add
+
+                    # Too low - made to avoid that atom is placed below slab
+                    if self.min_z is not None:
+                        tl = atoms_mutated.positions[i][2] < self.min_z
+                        if tl:
+                            continue
+                    else:
+                        tl = False
+                    
+                    # too close
+                    tc = atoms_too_close_two_sets(atoms_mutated[:i] + atoms_mutated[i+1:],
+                                         atoms_mutated[i], self.blmin)
+                    # too far
+                    if self.use_max_dist:
+                        if self.cluster:
+                            index_other = np.delete(np.arange(Nslab, Natoms), i)
+                        else:
+                            index_other = np.delete(np.arange(Natoms), i)
+                        blmin_array = np.array([self.blmin[num[i], num[j]] for j in index_other])
+                        tf = np.min(atoms_mutated.get_distances(i, index_other,mic=True)/blmin_array) > 1.7
+                    else:
+                        tf = False
+                    
+                    if not tc and not tf:
+                        break
+                if tc or tf or tl:
+                    atoms_mutated.positions[i] = posi_0
+        return atoms_mutated
+
+
+class RattleMutation_new(OffspringCreator):
+    """ 
+    A rattle mutation that that moves one atom and
+    and places it in the vicinity of one of the other
+    "top" atoms, randomly choosen.
+
+    Parameters:
+
+    blmin: Dictionary defining the minimum distance between atoms
+    after the rattle.
+    n_top: Number of atoms optimized by the GA.
+    rattle_prop: The probability with which each atom is rattled.
+    
+    """
+    def __init__(self, blmin, n_top, rattle_prop=0.4, min_z=None,
+                 verbose=False, descriptor='RattleMutation_new'):
+        OffspringCreator.__init__(self, verbose)
+        self.blmin = blmin
+        self.n_top = n_top
+        self.rattle_prop = rattle_prop
+        self.min_z = min_z
+        self.descriptor = descriptor
+        self.min_inputs = 1
+
+    def get_new_individual(self, parents, **kwargs):
+        f = parents[0]
+
+        indi = self.mutate(f)
+        if indi is None:
+            return indi, 'mutation: rattle'
+            
+        indi = self.initialize_individual(f, indi)
+        #indi.info['data']['parents'] = [f.info['confid']]  # remove for GP-search
+
+        return self.finalize_individual(indi), 'mutation: rattle'
+
+    def mutate(self, atoms):
+        num = atoms.numbers
+        atoms_mutated = atoms.copy()
+        Natoms = len(atoms)
+        Nslab = Natoms - self.n_top
+        for i in range(Nslab, Natoms):
+            if np.random.random() < self.rattle_prop:
+                posi_0 = np.copy(atoms_mutated.positions[i])
+                tc = False
+                tf = False
+
+                # Random order of atoms
+                #atom_indicies = np.random.permutation(np.arange(Nslab, Natoms))
+                atom_indicies = np.random.permutation(np.arange(Natoms))
+                for j in atom_indicies:
+                    if j==i:
+                        continue
+                    if self.min_z is not None:
+                        if atoms_mutated.positions[j][2] < self.min_z:
+                            continue
+                    posi_j = np.copy(atoms_mutated.positions[j])
+                    r_min = self.blmin[(num[i], num[j])]
+                    r_max = 1.7*r_min
+                    for k in range(20):
+                        atoms_mutated.positions[i] = posi_j
+                        # Then Rattle within a circle
+                        
+                        r = np.random.uniform(r_min**3, r_max**3)**(1/3)
+                        theta = np.random.uniform(low=0, high=2*np.pi)
+                        phi = np.random.uniform(low=0, high=np.pi)
+                        pos_add = r * np.array([np.cos(theta)*np.sin(phi),
+                                                np.sin(theta)*np.sin(phi),
+                                                np.cos(phi)])
+                        atoms_mutated.positions[i] += pos_add
+
+                        # Too low - made to avoid that atom is placed below slab
+                        if self.min_z is not None:
+                            tl = atoms_mutated.positions[i][2] < self.min_z
+                            if tl:
+                                continue
+                        else:
+                            tl = False
+                        
+                        # too close
+                        tc = atoms_too_close_two_sets(atoms_mutated[:i] + atoms_mutated[i+1:],
+                                                      atoms_mutated[i], self.blmin)
+                        # too far
+                        index_other = np.delete(np.arange(Natoms), i)
+                        blmin_array = np.array([self.blmin[num[i], num[j]] for j in index_other])
+                        tf = np.min(atoms_mutated.get_distances(i, index_other, mic=True)/blmin_array) > 1.7
+                        
+                        if not tc and not tf:
+                            break
+                    if tc or tf or tl:
+                        atoms_mutated.positions[i] = posi_0
+        return atoms_mutated
+    
+class PermutationMutation(OffspringCreator):
+    """Mutation that permutes a percentage of the atom types in the cluster.
+
+       Parameters:
+
+       n_top: Number of atoms optimized by the GA.
+       probability: The probability with which an atom is permuted.
+    """
+
+    def __init__(self, blmin, n_top, probability=0.33, verbose=False, descriptor='PermutationMutation'):
+        OffspringCreator.__init__(self, verbose)
+        self.blmin = blmin
+        self.n_top = n_top
+        self.probability = probability
+        self.descriptor = descriptor
+        self.min_inputs = 1
+
+    def get_new_individual(self, parents, **kwargs):
+        f = parents[0]
+
+        indi = self.mutate(f)
+        if indi is None:
+            return indi, 'mutation: permutation'
+            
+        indi = self.initialize_individual(f, indi)
+        # indi.info['data']['parents'] = [f.info['confid']]  # remove for GP-search
+
+        return self.finalize_individual(indi), 'mutation: permutation'
+
+    def mutate(self, atoms):
+        """ Does the actual mutation. """
+        a = atoms.copy()
+        Natoms = len(a)
+        Nslab = Natoms - self.n_top
+        pos = a.get_positions()[-self.n_top:]
+        num = a.numbers
+        num_top = num[-self.n_top:]
+        num_unique = list(set(num_top))
+        assert len(num_unique) > 1, 'Permutations with one atomic type is not valid'
+        m = int(ceil(float(self.n_top) * self.probability / 2.))
+        for _ in range(m):
+            swap_succesfull = False
+            for i in range(100):
+                # Find index pair of atoms to swap
+                i = j = 0
+                while num[i] == num[j]:
+                    i = np.random.randint(Nslab, Natoms)
+                    j = np.random.randint(Nslab, Natoms)
+
+                # Swap two atoms
+                pos_i = a.positions[i].copy()
+                pos_j = a.positions[j].copy()
+                a.positions[i] = pos_j
+                a.positions[j] = pos_i
+
+                # Check if atoms are too close
+                tc_i = atoms_too_close_two_sets(a[:i] + a[i+1:], a[i], self.blmin)
+                tc_j = atoms_too_close_two_sets(a[:j] + a[j+1:], a[j], self.blmin)
+                if tc_i or tc_j:
+                    # reset swap
+                    a.positions[i] = pos_i
+                    a.positions[j] = pos_j
+                    continue
+                else:
+                    swap_succesfull = True
+                    break
+            if not swap_succesfull:
+                # swap_back
+                a.positions[i] = pos_i
+                a.positions[j] = pos_j
+        return a
+
+
+class MirrorMutation(OffspringCreator):
+    """ A mirror mutation, as described in
+        TO BE PUBLISHED.
+        This mutation mirrors half of the cluster in a
+        randomly oriented cutting plane discarding the other half.
+
+        Parameters:
+        blmin: Dictionary defining the minimum allowed
+        distance between atoms.
+        n_top: Number of atoms the GA optimizes.
+        reflect: Defines if the mirrored half is also reflected
+        perpendicular to the mirroring plane.
+
+    """
+    def __init__(self, blmin, n_top, reflect=False, verbose=False):
+        OffspringCreator.__init__(self, verbose)
+        self.blmin = blmin
+        self.n_top = n_top
+        self.reflect = reflect
+        self.descriptor = 'MirrorMutation'
+        self.min_inputs = 1
+
+    def get_new_individual(self, parents):
+        f = parents[0]
+
+        indi = self.mutate(f)
+        if indi is None:
+            return indi, 'mutation: mirror'
+            
+        indi = self.initialize_individual(f, indi)
+        indi.info['data']['parents'] = [f.info['confid']]
+
+        return self.finalize_individual(indi), 'mutation: mirror'
+
+    def mutate(self, atoms):
+        """ Do the mutation of the atoms input. """
+
+        reflect = self.reflect
+        tc = True
+        slab = atoms[0:len(atoms) - self.n_top]
+        top = atoms[len(atoms) - self.n_top: len(atoms)]
+        num = top.numbers
+        unique_types = list(set(num))
+        nu = dict()
+        for u in unique_types:
+            nu[u] = sum(num == u)
+            
+        n_tries = 1000
+        counter = 0
+        changed = False
+
+        while tc and counter < n_tries:
+            counter += 1
+            cand = top.copy()
+            pos = cand.get_positions()
+
+            cm = np.average(top.get_positions(), axis=0)
+
+            # first select a randomly oriented cutting plane
+            theta = pi * np.random.random()
+            phi = 2. * pi * np.random.random()
+            n = (cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta))
+            n = np.array(n)
+
+            # Calculate all atoms signed distance to the cutting plane
+            D = []
+            for (i, p) in enumerate(pos):
+                d = np.dot(p - cm, n)
+                D.append((i, d))
+
+            # Sort the atoms by their signed distance
+            D.sort(key=lambda x: x[1])
+            nu_taken = dict()
+
+            # Select half of the atoms needed for