Commit d83d12b9 authored by Malthe Kjær Bisbo's avatar Malthe Kjær Bisbo
Browse files

dftb-calculator workaround, info-dict, initial structures now saved in all instances

parent b55781ba
......@@ -240,24 +240,24 @@ class PermutationMutation(CandidateGenerator):
indices_to_permute = [np.random.randint(Nslab,Natoms)]
# Perform permutations in sequence.
for i_permute in indices_to_permute:
for i in indices_to_permute:
for _ in range(100):
j_permute = np.random.randint(Nslab,Natoms)
while num[i_permute] == num[j_permute]:
j_permute = np.random.randint(Nslab,Natoms)
j = np.random.randint(Nslab,Natoms)
while num[i] == num[j]:
j = np.random.randint(Nslab,Natoms)
# Permute
pos_i = np.copy(a.positions[i_permute])
pos_j = np.copy(a.positions[j_permute])
a.positions[i_permute] = pos_j
a.positions[j_permute] = pos_i
pos_i = np.copy(a.positions[i])
pos_j = np.copy(a.positions[j])
a.positions[i] = pos_j
a.positions[j] = pos_i
# Check if rattle was valid
valid_bondlengths = self.check_valid_bondlengths(a, indices=[i,j])
if not valid_bondlengths:
a.positions[i_permute] = pos_i
a.positions[j_permute] = pos_j
a.positions[i] = pos_i
a.positions[j] = pos_j
else:
break
if valid_bondlengths:
......
"""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, oplist, probabilities_extra=None, oplist_extra=[], debugger=None):
assert len(probabilities) == len(oplist)
self.oplist = oplist
self.probabilities_extra = probabilities_extra
self.oplist_extra = oplist_extra
self.rho = np.cumsum(probabilities)
self.rho_extra = np.cumsum(probabilities_extra)
self.oplist_total = oplist + oplist_extra
self.debugger = debugger
def __get_index__(self, extra=False):
if extra:
v = np.random.random()
if v < self.rho_extra[-1]:
for i in range(len(self.rho_extra)):
if self.rho_extra[i] > v:
return i
else:
return None
else:
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, debugger=debugger)
if a is not None and self.probabilities_extra is not None:
to_use_extra = self.__get_index__(extra=True)
if to_use_extra is not None:
a, _ = self.oplist_extra[to_use_extra].get_new_individual([a,a])
return a, mut
else:
return a, mut
def get_operator(self):
"""Choose operator and return it."""
to_use = self.__get_index__()
return self.oplist[to_use]
""" 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 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]