Skip to main content
`Ctrl`+`K`

#

# Source code for ase.cluster.wulff

    # fmt: off

    import numpy as np

    delta = 1e-10

    def wulff_construction(symbol, surfaces, energies, size, structure,
                           rounding='closest', latticeconstant=None,
                           debug=False, maxiter=100):
        """Create a cluster using the Wulff construction.

        A cluster is created with approximately the number of atoms
        specified, following the Wulff construction, i.e. minimizing the
        surface energy of the cluster.

        Parameters
        ----------
        symbol : str or int
            The chemical symbol (or atomic number) of the desired element.

        surfaces : list
            A list of surfaces. Each surface is an (h, k, l) tuple or list of
            integers.

        energies : list
            A list of surface energies for the surfaces.

        size : int
            The desired number of atoms.

        structure : {'fcc', bcc', 'sc'}
            The desired crystal structure.

        rounding : {'closest', 'above', 'below'}
            Specifies what should be done if no Wulff construction corresponds
            to exactly the requested number of atoms. 'above', 'below', and
            'closest' mean that the nearest cluster above or below - or the
            closest one - is created instead.

        latticeconstant : float (optional)
            The lattice constant. If not given, extracted from `ase.data`.

        debug : bool, default False
            If True, information about the iteration towards the right cluster
            size is printed.
        """

        if debug:
            print('Wulff: Aiming for cluster with %i atoms (%s)' %
                  (size, rounding))

            if rounding not in ['above', 'below', 'closest']:
                raise ValueError(f'Invalid rounding: {rounding}')

        # Interpret structure, if it is a string.
        if isinstance(structure, str):
            if structure == 'fcc':
                from ase.cluster.cubic import FaceCenteredCubic as structure
            elif structure == 'bcc':
                from ase.cluster.cubic import BodyCenteredCubic as structure
            elif structure == 'sc':
                from ase.cluster.cubic import SimpleCubic as structure
            elif structure == 'hcp':
                from ase.cluster.hexagonal import HexagonalClosedPacked as structure
            elif structure == 'graphite':
                from ase.cluster.hexagonal import Graphite as structure
            else:
                error = f'Crystal structure {structure} is not supported.'
                raise NotImplementedError(error)

        # Check number of surfaces
        nsurf = len(surfaces)
        if len(energies) != nsurf:
            raise ValueError('The energies array should contain %d values.'
                             % (nsurf,))

        # Copy energies array so it is safe to modify it
        energies = np.array(energies)

        # We should check that for each direction, the surface energy plus
        # the energy in the opposite direction is positive.  But this is
        # very difficult in the general case!

        # Before starting, make a fake cluster just to extract the
        # interlayer distances in the relevant directions, and use these
        # to "renormalize" the surface energies such that they can be used
        # to convert to number of layers instead of to distances.
        atoms = structure(symbol, surfaces, 5 * np.ones(len(surfaces), int),
                          latticeconstant=latticeconstant)
        for i, s in enumerate(surfaces):
            d = atoms.get_layer_distance(s)
            energies[i] /= d

        # First guess a size that is not too large.
        wanted_size = size ** (1.0 / 3.0)
        max_e = max(energies)
        factor = wanted_size / max_e
        atoms, layers = make_atoms(symbol, surfaces, energies, factor, structure,
                                   latticeconstant)
        if len(atoms) == 0:
            # Probably the cluster is very flat
            if debug:
                print('First try made an empty cluster, trying again.')
            factor = 1 / energies.min()
            atoms, layers = make_atoms(symbol, surfaces, energies, factor,
                                       structure, latticeconstant)
            if len(atoms) == 0:
                raise RuntimeError('Failed to create a finite cluster.')

        # Second guess: scale to get closer.
        old_factor = factor
        old_layers = layers
        old_atoms = atoms
        factor *= (size / len(atoms))**(1.0 / 3.0)
        atoms, layers = make_atoms(symbol, surfaces, energies, factor,
                                   structure, latticeconstant)
        if len(atoms) == 0:
            print('Second guess gave an empty cluster, discarding it.')
            atoms = old_atoms
            factor = old_factor
            layers = old_layers
        else:
            del old_atoms

        # Find if the cluster is too small or too large (both means perfect!)
        below = above = None
        if len(atoms) <= size:
            below = atoms
        if len(atoms) >= size:
            above = atoms

        # Now iterate towards the right cluster
        iter = 0
        while (below is None or above is None):
            if len(atoms) < size:
                # Find a larger cluster
                if debug:
                    print('Making a larger cluster.')
                factor = ((layers + 0.5 + delta) / energies).min()
                atoms, new_layers = make_atoms(symbol, surfaces, energies, factor,
                                               structure, latticeconstant)
                assert (new_layers - layers).max() == 1
                assert (new_layers - layers).min() >= 0
                layers = new_layers
            else:
                # Find a smaller cluster
                if debug:
                    print('Making a smaller cluster.')
                factor = ((layers - 0.5 - delta) / energies).max()
                atoms, new_layers = make_atoms(symbol, surfaces, energies, factor,
                                               structure, latticeconstant)
                assert (new_layers - layers).max() <= 0
                assert (new_layers - layers).min() == -1
                layers = new_layers
            if len(atoms) <= size:
                below = atoms
            if len(atoms) >= size:
                above = atoms
            iter += 1
            if iter == maxiter:
                raise RuntimeError('Runaway iteration.')
        if rounding == 'below':
            if debug:
                print('Choosing smaller cluster with %i atoms' % len(below))
            return below
        elif rounding == 'above':
            if debug:
                print('Choosing larger cluster with %i atoms' % len(above))
            return above
        else:
            assert rounding == 'closest'
            if (len(above) - size) < (size - len(below)):
                atoms = above
            else:
                atoms = below
            if debug:
                print('Choosing closest cluster with %i atoms' % len(atoms))
            return atoms

    def make_atoms(symbol, surfaces, energies, factor, structure, latticeconstant):
        layers1 = factor * np.array(energies)
        layers = np.round(layers1).astype(int)
        atoms = structure(symbol, surfaces, layers,
                          latticeconstant=latticeconstant)
        return (atoms, layers)