Source code for ase.ga.standardmutations

"""A collection of mutations that can be used."""
import numpy as np
from math import cos, sin, pi
from ase.calculators.lammpslib import convert_cell
from ase.ga.utilities import (atoms_too_close,
                              atoms_too_close_two_sets,
                              gather_atoms_by_tag,
                              get_rotation_matrix)
from ase.ga.offspring_creator import OffspringCreator, CombinationMutation
from ase import Atoms


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.

    test_dist_to_slab: whether to also make sure that the distances
        between the atoms and the slab satisfy the blmin.

    use_tags: if True, the atomic tags will be used to preserve
        molecular identity. Same-tag atoms will then be
        displaced collectively, so that the internal
        geometry is preserved.

    rng: Random number generator
        By default numpy.random.
    """
    def __init__(self, blmin, n_top, rattle_strength=0.8,
                 rattle_prop=0.4, test_dist_to_slab=True, use_tags=False,
                 verbose=False, rng=np.random):
        OffspringCreator.__init__(self, verbose, rng=rng)
        self.blmin = blmin
        self.n_top = n_top
        self.rattle_strength = rattle_strength
        self.rattle_prop = rattle_prop
        self.test_dist_to_slab = test_dist_to_slab
        self.use_tags = use_tags

        self.descriptor = 'RattleMutation'
        self.min_inputs = 1

    def get_new_individual(self, parents):
        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']]

        return self.finalize_individual(indi), 'mutation: rattle'

    def mutate(self, atoms):
        """Does the actual mutation."""
        N = len(atoms) if self.n_top is None else self.n_top
        slab = atoms[:len(atoms) - N]
        atoms = atoms[-N:]
        tags = atoms.get_tags() if self.use_tags else np.arange(N)
        pos_ref = atoms.get_positions()
        num = atoms.get_atomic_numbers()
        cell = atoms.get_cell()
        pbc = atoms.get_pbc()
        st = 2. * self.rattle_strength

        count = 0
        maxcount = 1000
        too_close = True
        while too_close and count < maxcount:
            count += 1
            pos = pos_ref.copy()
            ok = False
            for tag in np.unique(tags):
                select = np.where(tags == tag)
                if self.rng.rand() < self.rattle_prop:
                    ok = True
                    r = self.rng.rand(3)
                    pos[select] += st * (r - 0.5)

            if not ok:
                # Nothing got rattled
                continue

            top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
            too_close = atoms_too_close(
                top, self.blmin, use_tags=self.use_tags)
            if not too_close and self.test_dist_to_slab:
                too_close = atoms_too_close_two_sets(top, slab, self.blmin)

        if count == maxcount:
            return None

        mutant = slab + top
        return mutant


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.

    test_dist_to_slab: whether to also make sure that the distances
        between the atoms and the slab satisfy the blmin.

    use_tags: if True, the atomic tags will be used to preserve
        molecular identity. Permutations will then happen
        at the molecular level, i.e. swapping the center-of-
        positions of two moieties while preserving their
        internal geometries.

    blmin: Dictionary defining the minimum distance between atoms
        after the permutation. If equal to None (the default),
        no such check is performed.

    rng: Random number generator
        By default numpy.random.
    """
    def __init__(self, n_top, probability=0.33, test_dist_to_slab=True,
                 use_tags=False, blmin=None, rng=np.random, verbose=False):
        OffspringCreator.__init__(self, verbose, rng=rng)
        self.n_top = n_top
        self.probability = probability
        self.test_dist_to_slab = test_dist_to_slab
        self.use_tags = use_tags
        self.blmin = blmin

        self.descriptor = 'PermutationMutation'
        self.min_inputs = 1

    def get_new_individual(self, parents):
        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']]

        return self.finalize_individual(indi), 'mutation: permutation'

    def mutate(self, atoms):
        """Does the actual mutation."""
        N = len(atoms) if self.n_top is None else self.n_top
        slab = atoms[:len(atoms) - N]
        atoms = atoms[-N:]
        if self.use_tags:
            gather_atoms_by_tag(atoms)
        tags = atoms.get_tags() if self.use_tags else np.arange(N)
        pos_ref = atoms.get_positions()
        num = atoms.get_atomic_numbers()
        cell = atoms.get_cell()
        pbc = atoms.get_pbc()
        symbols = atoms.get_chemical_symbols()

        unique_tags = np.unique(tags)
        n = len(unique_tags)
        swaps = int(np.ceil(n * self.probability / 2.))

        sym = []
        for tag in unique_tags:
            indices = np.where(tags == tag)[0]
            s = ''.join([symbols[j] for j in indices])
            sym.append(s)

        assert len(np.unique(sym)) > 1, \
            'Permutations with one atom (or molecule) type is not valid'

        count = 0
        maxcount = 1000
        too_close = True
        while too_close and count < maxcount:
            count += 1
            pos = pos_ref.copy()
            for _ in range(swaps):
                i = j = 0
                while sym[i] == sym[j]:
                    i = self.rng.randint(0, high=n)
                    j = self.rng.randint(0, high=n)
                ind1 = np.where(tags == i)
                ind2 = np.where(tags == j)
                cop1 = np.mean(pos[ind1], axis=0)
                cop2 = np.mean(pos[ind2], axis=0)
                pos[ind1] += cop2 - cop1
                pos[ind2] += cop1 - cop2

            top = Atoms(num, positions=pos, cell=cell, pbc=pbc, tags=tags)
            if self.blmin is None:
                too_close = False
            else:
                too_close = atoms_too_close(
                    top, self.blmin, use_tags=self.use_tags)
                if not too_close and self.test_dist_to_slab:
                    too_close = atoms_too_close_two_sets(top, slab, self.blmin)

        if count == maxcount:
            return None

        mutant = slab + top
        return mutant


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.

    rng: Random number generator
        By default numpy.random.
    """
    def __init__(self, blmin, n_top, reflect=False, rng=np.random,
                 verbose=False):
        OffspringCreator.__init__(self, verbose, rng=rng)
        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 * self.rng.rand()
            phi = 2. * pi * self.rng.rand()
            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 a full cluster
            p_use = []
            n_use = []
            for (i, d) in D:
                if num[i] not in nu_taken.keys():
                    nu_taken[num[i]] = 0
                if nu_taken[num[i]] < nu[num[i]] / 2.:
                    p_use.append(pos[i])
                    n_use.append(num[i])
                    nu_taken[num[i]] += 1

            # calculate the mirrored position and add these.
            pn = []
            for p in p_use:
                pt = p - 2. * np.dot(p - cm, n) * n
                if reflect:
                    pt = -pt + 2 * cm + 2 * n * np.dot(pt - cm, n)
                pn.append(pt)

            n_use.extend(n_use)
            p_use.extend(pn)

            # In the case of an uneven number of
            # atoms we need to add one extra
            for n in nu.keys():
                if nu[n] % 2 == 0:
                    continue
                while sum(n_use == n) > nu[n]:
                    for i in range(int(len(n_use) / 2), len(n_use)):
                        if n_use[i] == n:
                            del p_use[i]
                            del n_use[i]
                            break
                assert sum(n_use == n) == nu[n]

            # Make sure we have the correct number of atoms
            # and rearrange the atoms so they are in the right order
            for i in range(len(n_use)):
                if num[i] == n_use[i]:
                    continue
                for j in range(i + 1, len(n_use)):
                    if n_use[j] == num[i]:
                        tn = n_use[i]
                        tp = p_use[i]
                        n_use[i] = n_use[j]
                        p_use[i] = p_use[j]
                        p_use[j] = tp
                        n_use[j] = tn

            # Finally we check that nothing is too close in the end product.
            cand = Atoms(num, p_use, cell=slab.get_cell(), pbc=slab.get_pbc())

            tc = atoms_too_close(cand, self.blmin)
            if tc:
                continue
            tc = atoms_too_close_two_sets(slab, cand, self.blmin)

            if not changed and counter > n_tries // 2:
                reflect = not reflect
                changed = True

            tot = slab + cand

        if counter == n_tries:
            return None
        return tot


[docs]class StrainMutation(OffspringCreator): """ Mutates a candidate by applying a randomly generated strain. For more information, see also: * `Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720`__ __ https://doi.org/10.1016/j.cpc.2006.07.020 * `Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387`__ __ https://doi.org/10.1016/j.cpc.2010.07.048 After initialization of the mutation, a scaling volume (to which each mutated structure is scaled before checking the constraints) is typically generated from the population, which is then also occasionally updated in the course of the GA run. Parameters: blmin: dict The closest allowed interatomic distances on the form: {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers. cellbounds: ase.ga.utilities.CellBounds instance Describes limits on the cell shape, see :class:`~ase.ga.utilities.CellBounds`. stddev: float Standard deviation used in the generation of the strain matrix elements. number_of_variable_cell_vectors: int (default 3) The number of variable cell vectors (1, 2 or 3). To keep things simple, it is the 'first' vectors which will be treated as variable, i.e. the 'a' vector in the univariate case, the 'a' and 'b' vectors in the bivariate case, etc. use_tags: boolean Whether to use the atomic tags to preserve molecular identity. rng: Random number generator By default numpy.random. """ def __init__(self, blmin, cellbounds=None, stddev=0.7, number_of_variable_cell_vectors=3, use_tags=False, rng=np.random, verbose=False): OffspringCreator.__init__(self, verbose, rng=rng) self.blmin = blmin self.cellbounds = cellbounds self.stddev = stddev self.number_of_variable_cell_vectors = number_of_variable_cell_vectors self.use_tags = use_tags self.scaling_volume = None self.descriptor = 'StrainMutation' self.min_inputs = 1 def update_scaling_volume(self, population, w_adapt=0.5, n_adapt=0): """Function to initialize or update the scaling volume in a GA run. w_adapt: weight of the new vs the old scaling volume n_adapt: number of best candidates in the population that are used to calculate the new scaling volume """ if not n_adapt: # if not set, take best 20% of the population n_adapt = int(np.ceil(0.2 * len(population))) v_new = np.mean([a.get_volume() for a in population[:n_adapt]]) if not self.scaling_volume: self.scaling_volume = v_new else: volumes = [self.scaling_volume, v_new] weights = [1 - w_adapt, w_adapt] self.scaling_volume = np.average(volumes, weights=weights) def get_new_individual(self, parents): f = parents[0] indi = self.mutate(f) if indi is None: return indi, 'mutation: strain' indi = self.initialize_individual(f, indi) indi.info['data']['parents'] = [f.info['confid']] return self.finalize_individual(indi), 'mutation: strain' def mutate(self, atoms): """ Does the actual mutation. """ cell_ref = atoms.get_cell() pos_ref = atoms.get_positions() vol_ref = atoms.get_volume() if self.use_tags: tags = atoms.get_tags() gather_atoms_by_tag(atoms) pos = atoms.get_positions() mutant = atoms.copy() count = 0 too_close = True maxcount = 1000 while too_close and count < maxcount: count += 1 # generating the strain matrix: strain = np.identity(3) for i in range(self.number_of_variable_cell_vectors): for j in range(i + 1): r = self.rng.normal(loc=0., scale=self.stddev) if i == j: strain[i, j] += r else: epsilon = 0.5 * r strain[i, j] += epsilon strain[j, i] += epsilon # applying the strain: cell_new = np.dot(strain, cell_ref) # convert to lower triangular form cell_new = convert_cell(cell_new)[0].T # volume scaling: if self.number_of_variable_cell_vectors > 0: volume = abs(np.linalg.det(cell_new)) if self.scaling_volume is None: # The scaling_volume has not been set (yet), # so we give it the same volume as the parent scaling = vol_ref / volume else: scaling = self.scaling_volume / volume scaling **= 1. / self.number_of_variable_cell_vectors cell_new[:self.number_of_variable_cell_vectors] *= scaling # check cell dimensions: if not self.cellbounds.is_within_bounds(cell_new): continue # ensure non-variable cell vectors are indeed unchanged for i in range(self.number_of_variable_cell_vectors, 3): assert np.allclose(cell_new[i], cell_ref[i]) # apply the new unit cell and scale # the atomic positions accordingly mutant.set_cell(cell_ref, scale_atoms=False) if self.use_tags: transfo = np.linalg.solve(cell_ref, cell_new) for tag in np.unique(tags): select = np.where(tags == tag) cop = np.mean(pos[select], axis=0) disp = np.dot(cop, transfo) - cop mutant.positions[select] += disp else: mutant.set_positions(pos_ref) mutant.set_cell(cell_new, scale_atoms=not self.use_tags) mutant.wrap() # check the interatomic distances too_close = atoms_too_close(mutant, self.blmin, use_tags=self.use_tags) if count == maxcount: mutant = None return mutant
[docs]class PermuStrainMutation(CombinationMutation): """Combination of PermutationMutation and StrainMutation. For more information, see also: * `Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387`__ __ https://doi.org/10.1016/j.cpc.2010.07.048 Parameters: permutationmutation: OffspringCreator instance A mutation that permutes atom types. strainmutation: OffspringCreator instance A mutation that mutates by straining. """ def __init__(self, permutationmutation, strainmutation, verbose=False): super(PermuStrainMutation, self).__init__(permutationmutation, strainmutation, verbose=verbose) self.descriptor = 'permustrain'
[docs]class RotationalMutation(OffspringCreator): """Mutates a candidate by applying random rotations to multi-atom moieties in the structure (atoms with the same tag are considered part of one such moiety). Only performs whole-molecule rotations, no internal rotations. For more information, see also: * `Zhu Q., Oganov A.R., Glass C.W., Stokes H.T, Acta Cryst. (2012), B68, 215-226.`__ __ https://dx.doi.org/10.1107/S0108768112017466 Parameters: blmin: dict The closest allowed interatomic distances on the form: {(Z, Z*): dist, ...}, where Z and Z* are atomic numbers. n_top: int or None The number of atoms to optimize (None = include all). fraction: float Fraction of the moieties to be rotated. tags: None or list of integers Specifies, respectively, whether all moieties or only those with matching tags are eligible for rotation. min_angle: float Minimal angle (in radians) for each rotation; should lie in the interval [0, pi]. test_dist_to_slab: boolean Whether also the distances to the slab should be checked to satisfy the blmin. rng: Random number generator By default numpy.random. """ def __init__(self, blmin, n_top=None, fraction=0.33, tags=None, min_angle=1.57, test_dist_to_slab=True, rng=np.random, verbose=False): OffspringCreator.__init__(self, verbose, rng=rng) self.blmin = blmin self.n_top = n_top self.fraction = fraction self.tags = tags self.min_angle = min_angle self.test_dist_to_slab = test_dist_to_slab self.descriptor = 'RotationalMutation' self.min_inputs = 1 def get_new_individual(self, parents): f = parents[0] indi = self.mutate(f) if indi is None: return indi, 'mutation: rotational' indi = self.initialize_individual(f, indi) indi.info['data']['parents'] = [f.info['confid']] return self.finalize_individual(indi), 'mutation: rotational' def mutate(self, atoms): """Does the actual mutation.""" N = len(atoms) if self.n_top is None else self.n_top slab = atoms[:len(atoms) - N] atoms = atoms[-N:] mutant = atoms.copy() gather_atoms_by_tag(mutant) pos = mutant.get_positions() tags = mutant.get_tags() eligible_tags = tags if self.tags is None else self.tags indices = {} for tag in np.unique(tags): hits = np.where(tags == tag)[0] if len(hits) > 1 and tag in eligible_tags: indices[tag] = hits n_rot = int(np.ceil(len(indices) * self.fraction)) chosen_tags = self.rng.choice(list(indices.keys()), size=n_rot, replace=False) too_close = True count = 0 maxcount = 10000 while too_close and count < maxcount: newpos = np.copy(pos) for tag in chosen_tags: p = np.copy(newpos[indices[tag]]) cop = np.mean(p, axis=0) if len(p) == 2: line = (p[1] - p[0]) / np.linalg.norm(p[1] - p[0]) while True: axis = self.rng.rand(3) axis /= np.linalg.norm(axis) a = np.arccos(np.dot(axis, line)) if np.pi / 4 < a < np.pi * 3 / 4: break else: axis = self.rng.rand(3) axis /= np.linalg.norm(axis) angle = self.min_angle angle += 2 * (np.pi - self.min_angle) * self.rng.rand() m = get_rotation_matrix(axis, angle) newpos[indices[tag]] = np.dot(m, (p - cop).T).T + cop mutant.set_positions(newpos) mutant.wrap() too_close = atoms_too_close(mutant, self.blmin, use_tags=True) count += 1 if not too_close and self.test_dist_to_slab: too_close = atoms_too_close_two_sets(slab, mutant, self.blmin) if count == maxcount: mutant = None else: mutant = slab + mutant return mutant
[docs]class RattleRotationalMutation(CombinationMutation): """Combination of RattleMutation and RotationalMutation. Parameters: rattlemutation: OffspringCreator instance A mutation that rattles atoms. rotationalmutation: OffspringCreator instance A mutation that rotates moieties. """ def __init__(self, rattlemutation, rotationalmutation, verbose=False): super(RattleRotationalMutation, self).__init__(rattlemutation, rotationalmutation, verbose=verbose) self.descriptor = 'rattlerotational'