"""Implementation of the cut-and-splice paring operator."""
import numpy as np
from ase import Atoms
from ase.ga.offspring_creator import OffspringCreator
from ase.ga.utilities import (
atoms_too_close,
atoms_too_close_two_sets,
gather_atoms_by_tag,
)
from ase.geometry import find_mic
class Positions:
"""Helper object to simplify the pairing process.
Parameters:
scaled_positions: (Nx3) array
Positions in scaled coordinates
cop: (1x3) array
Center-of-positions (also in scaled coordinates)
symbols: str
String with the atomic symbols
distance: float
Signed distance to the cutting plane
origin: int (0 or 1)
Determines at which side of the plane the position should be.
"""
def __init__(self, scaled_positions, cop, symbols, distance, origin):
self.scaled_positions = scaled_positions
self.cop = cop
self.symbols = symbols
self.distance = distance
self.origin = origin
def to_use(self):
"""Tells whether this position is at the right side."""
if self.distance > 0. and self.origin == 0:
return True
elif self.distance < 0. and self.origin == 1:
return True
else:
return False
[docs]
class CutAndSplicePairing(OffspringCreator):
"""The Cut and Splice operator by Deaven and Ho.
Creates offspring from two parent structures using
a randomly generated cutting plane.
The parents may have different unit cells, in which
case the offspring unit cell will be a random combination
of the parent cells.
The basic implementation (for fixed unit cells) is
described in:
:doi:`L.B. Vilhelmsen and B. Hammer, PRL, 108, 126101 (2012)
<10.1103/PhysRevLett.108.126101`>
The extension to variable unit cells is similar to:
* :doi:`Glass, Oganov, Hansen, Comp. Phys. Comm. 175 (2006) 713-720
<10.1016/j.cpc.2006.07.020>`
* :doi:`Lonie, Zurek, Comp. Phys. Comm. 182 (2011) 372-387
<10.1016/j.cpc.2010.07.048>`
The operator can furthermore preserve molecular identity
if desired (see the *use_tags* kwarg). Atoms with the same
tag will then be considered as belonging to the same molecule,
and their internal geometry will not be changed by the operator.
If use_tags is enabled, the operator will also conserve the
number of molecules of each kind (in addition to conserving
the overall stoichiometry). Currently, molecules are considered
to be of the same kind if their chemical symbol strings are
identical. In rare cases where this may not be sufficient
(i.e. when desiring to keep the same ratio of isomers), the
different isomers can be differentiated by giving them different
elemental orderings (e.g. 'XY2' and 'Y2X').
Parameters:
slab: Atoms object
Specifies the cell vectors and periodic boundary conditions
to be applied to the randomly generated structures.
Any included atoms (e.g. representing an underlying slab)
are copied to these new structures.
n_top: int
The number of atoms to optimize
blmin: dict
Dictionary with minimal interatomic distances.
Note: when preserving molecular identity (see use_tags),
the blmin dict will (naturally) only be applied
to intermolecular distances (not the intramolecular
ones).
number_of_variable_cell_vectors: int (default 0)
The number of variable cell vectors (0, 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.
p1: float or int between 0 and 1
Probability that a parent is shifted over a random
distance along the normal of the cutting plane
(only operative if number_of_variable_cell_vectors > 0).
p2: float or int between 0 and 1
Same as p1, but for shifting along the directions
in the cutting plane (only operative if
number_of_variable_cell_vectors > 0).
minfrac: float between 0 and 1, or None (default)
Minimal fraction of atoms a parent must contribute
to the child. If None, each parent must contribute
at least one atom.
cellbounds: ase.ga.utilities.CellBounds instance
Describing limits on the cell shape, see
:class:`~ase.ga.utilities.CellBounds`.
Note that it only make sense to impose conditions
regarding cell vectors which have been marked as
variable (see number_of_variable_cell_vectors).
use_tags: bool
Whether to use the atomic tags to preserve
molecular identity.
test_dist_to_slab: bool (default True)
Whether to make sure that the distances between
the atoms and the slab satisfy the blmin.
rng: Random number generator
By default numpy.random.
"""
def __init__(self, slab, n_top, blmin, number_of_variable_cell_vectors=0,
p1=1, p2=0.05, minfrac=None, cellbounds=None,
test_dist_to_slab=True, use_tags=False, rng=np.random,
verbose=False):
OffspringCreator.__init__(self, verbose, rng=rng)
self.slab = slab
self.n_top = n_top
self.blmin = blmin
assert number_of_variable_cell_vectors in range(4)
self.number_of_variable_cell_vectors = number_of_variable_cell_vectors
self.p1 = p1
self.p2 = p2
self.minfrac = minfrac
self.cellbounds = cellbounds
self.test_dist_to_slab = test_dist_to_slab
self.use_tags = use_tags
self.scaling_volume = None
self.descriptor = 'CutAndSplicePairing'
self.min_inputs = 2
def update_scaling_volume(self, population, w_adapt=0.5, n_adapt=0):
"""Updates the scaling volume that is used in the pairing
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:
# 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):
"""The method called by the user that
returns the paired structure."""
f, m = parents
indi = self.cross(f, m)
desc = 'pairing: {} {}'.format(f.info['confid'],
m.info['confid'])
# It is ok for an operator to return None
# It means that it could not make a legal offspring
# within a reasonable amount of time
if indi is None:
return indi, desc
indi = self.initialize_individual(f, indi)
indi.info['data']['parents'] = [f.info['confid'],
m.info['confid']]
return self.finalize_individual(indi), desc
def cross(self, a1, a2):
"""Crosses the two atoms objects and returns one"""
if len(a1) != len(self.slab) + self.n_top:
raise ValueError('Wrong size of structure to optimize')
if len(a1) != len(a2):
raise ValueError('The two structures do not have the same length')
N = self.n_top
# Only consider the atoms to optimize
a1 = a1[len(a1) - N: len(a1)]
a2 = a2[len(a2) - N: len(a2)]
if not np.array_equal(a1.numbers, a2.numbers):
err = 'Trying to pair two structures with different stoichiometry'
raise ValueError(err)
if self.use_tags and not np.array_equal(a1.get_tags(), a2.get_tags()):
err = 'Trying to pair two structures with different tags'
raise ValueError(err)
cell1 = a1.get_cell()
cell2 = a2.get_cell()
for i in range(self.number_of_variable_cell_vectors, 3):
err = 'Unit cells are supposed to be identical in direction %d'
assert np.allclose(cell1[i], cell2[i]), (err % i, cell1, cell2)
invalid = True
counter = 0
maxcount = 1000
a1_copy = a1.copy()
a2_copy = a2.copy()
# Run until a valid pairing is made or maxcount pairings are tested.
while invalid and counter < maxcount:
counter += 1
newcell = self.generate_unit_cell(cell1, cell2)
if newcell is None:
# No valid unit cell could be generated.
# This strongly suggests that it is near-impossible
# to generate one from these parent cells and it is
# better to abort now.
break
# Choose direction of cutting plane normal
if self.number_of_variable_cell_vectors == 0:
# Will be generated entirely at random
theta = np.pi * self.rng.random()
phi = 2. * np.pi * self.rng.random()
cut_n = np.array([np.cos(phi) * np.sin(theta),
np.sin(phi) * np.sin(theta), np.cos(theta)])
else:
# Pick one of the 'variable' cell vectors
cut_n = self.rng.choice(self.number_of_variable_cell_vectors)
# Randomly translate parent structures
for a_copy, a in zip([a1_copy, a2_copy], [a1, a2]):
a_copy.set_positions(a.get_positions())
cell = a_copy.get_cell()
for i in range(self.number_of_variable_cell_vectors):
r = self.rng.random()
cond1 = i == cut_n and r < self.p1
cond2 = i != cut_n and r < self.p2
if cond1 or cond2:
a_copy.positions += self.rng.random() * cell[i]
if self.use_tags:
# For correct determination of the center-
# of-position of the multi-atom blocks,
# we need to group their constituent atoms
# together
gather_atoms_by_tag(a_copy)
else:
a_copy.wrap()
# Generate the cutting point in scaled coordinates
cosp1 = np.average(a1_copy.get_scaled_positions(), axis=0)
cosp2 = np.average(a2_copy.get_scaled_positions(), axis=0)
cut_p = np.zeros((1, 3))
for i in range(3):
if i < self.number_of_variable_cell_vectors:
cut_p[0, i] = self.rng.random()
else:
cut_p[0, i] = 0.5 * (cosp1[i] + cosp2[i])
# Perform the pairing:
child = self._get_pairing(a1_copy, a2_copy, cut_p, cut_n, newcell)
if child is None:
continue
# Verify whether the atoms are too close or not:
if atoms_too_close(child, self.blmin, use_tags=self.use_tags):
continue
if self.test_dist_to_slab and len(self.slab) > 0:
if atoms_too_close_two_sets(self.slab, child, self.blmin):
continue
# Passed all the tests
child = self.slab + child
child.set_cell(newcell, scale_atoms=False)
child.wrap()
return child
return None
def generate_unit_cell(self, cell1, cell2, maxcount=10000):
"""Generates a new unit cell by a random linear combination
of the parent cells. The new cell must satisfy the
self.cellbounds constraints. Returns None if no such cell
was generated within a given number of attempts.
Parameters:
maxcount: int
The maximal number of attempts.
"""
# First calculate the scaling volume
if not self.scaling_volume:
v1 = np.abs(np.linalg.det(cell1))
v2 = np.abs(np.linalg.det(cell2))
r = self.rng.random()
v_ref = r * v1 + (1 - r) * v2
else:
v_ref = self.scaling_volume
# Now the cell vectors
if self.number_of_variable_cell_vectors == 0:
assert np.allclose(cell1, cell2), 'Parent cells are not the same'
newcell = np.copy(cell1)
else:
count = 0
while count < maxcount:
r = self.rng.random()
newcell = r * cell1 + (1 - r) * cell2
vol = abs(np.linalg.det(newcell))
scaling = v_ref / vol
scaling **= 1. / self.number_of_variable_cell_vectors
newcell[:self.number_of_variable_cell_vectors] *= scaling
found = True
if self.cellbounds is not None:
found = self.cellbounds.is_within_bounds(newcell)
if found:
break
count += 1
else:
# Did not find acceptable cell
newcell = None
return newcell
def _get_pairing(self, a1, a2, cutting_point, cutting_normal, cell):
"""Creates a child from two parents using the given cut.
Returns None if the generated structure does not contain
a large enough fraction of each parent (see self.minfrac).
Does not check whether atoms are too close.
Assumes the 'slab' parts have been removed from the parent
structures and that these have been checked for equal
lengths, stoichiometries, and tags (if self.use_tags).
Parameters:
cutting_normal: int or (1x3) array
cutting_point: (1x3) array
In fractional coordinates
cell: (3x3) array
The unit cell for the child structure
"""
symbols = a1.get_chemical_symbols()
tags = a1.get_tags() if self.use_tags else np.arange(len(a1))
# Generate list of all atoms / atom groups:
p1, p2, sym = [], [], []
for i in np.unique(tags):
indices = np.where(tags == i)[0]
s = ''.join([symbols[j] for j in indices])
sym.append(s)
for i, (a, p) in enumerate(zip([a1, a2], [p1, p2])):
c = a.get_cell()
cop = np.mean(a.positions[indices], axis=0)
cut_p = np.dot(cutting_point, c)
if isinstance(cutting_normal, int):
vecs = [c[j] for j in range(3) if j != cutting_normal]
cut_n = np.cross(vecs[0], vecs[1])
else:
cut_n = np.dot(cutting_normal, c)
d = np.dot(cop - cut_p, cut_n)
spos = a.get_scaled_positions()[indices]
scop = np.mean(spos, axis=0)
p.append(Positions(spos, scop, s, d, i))
all_points = p1 + p2
unique_sym = np.unique(sym)
types = {s: sym.count(s) for s in unique_sym}
# Sort these by chemical symbols:
all_points.sort(key=lambda x: x.symbols, reverse=True)
# For each atom type make the pairing
unique_sym.sort()
use_total = {}
for s in unique_sym:
used = []
not_used = []
# The list is looked trough in
# reverse order so atoms can be removed
# from the list along the way.
for i in reversed(range(len(all_points))):
# If there are no more atoms of this type
if all_points[i].symbols != s:
break
# Check if the atom should be included
if all_points[i].to_use():
used.append(all_points.pop(i))
else:
not_used.append(all_points.pop(i))
assert len(used) + len(not_used) == types[s] * 2
# While we have too few of the given atom type
while len(used) < types[s]:
index = self.rng.randint(len(not_used))
used.append(not_used.pop(index))
# While we have too many of the given atom type
while len(used) > types[s]:
# remove randomly:
index = self.rng.randint(len(used))
not_used.append(used.pop(index))
use_total[s] = used
n_tot = sum(len(ll) for ll in use_total.values())
assert n_tot == len(sym)
# check if the generated structure contains
# atoms from both parents:
count1, count2, N = 0, 0, len(a1)
for x in use_total.values():
count1 += sum(y.origin == 0 for y in x)
count2 += sum(y.origin == 1 for y in x)
nmin = 1 if self.minfrac is None else int(round(self.minfrac * N))
if count1 < nmin or count2 < nmin:
return None
# Construct the cartesian positions and reorder the atoms
# to follow the original order
newpos = []
pbc = a1.get_pbc()
for s in sym:
p = use_total[s].pop()
c = a1.get_cell() if p.origin == 0 else a2.get_cell()
pos = np.dot(p.scaled_positions, c)
cop = np.dot(p.cop, c)
vectors, _lengths = find_mic(pos - cop, c, pbc)
newcop = np.dot(p.cop, cell)
pos = newcop + vectors
for row in pos:
newpos.append(row)
newpos = np.reshape(newpos, (N, 3))
num = a1.get_atomic_numbers()
child = Atoms(numbers=num, positions=newpos, pbc=pbc, cell=cell,
tags=tags)
child.wrap()
return child