Skip to content

GitLab

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

rattleMutation, startGenerator, OperationsSelector, CandidateGenerator working... 

rattleMutation, startGenerator, OperationsSelector, CandidateGenerator working - but test a bit more
parent a42a26e9
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 539 additions and 565 deletions
bfgslinesearch_constrained.py 0 → 100644
View file @ 5c688cd1
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
from utils import check_valid_bondlengths
# 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_constrained(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,
blmin=None, blmax=None, dmax_cov=3.0, 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.blmin = blmin
self.blmax = blmax
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 constraints
terminate_due_to_constraints = self.check_constraints()
if terminate_due_to_constraints:
return True
return (forces**2).sum(axis=1).max() < self.fmax**2
def check_constraints(self):
# Check if stop due to large displacement
valid_displace = self.check_displacement()
# Check if stop due to zlim
valid_z = self.check_zcoord()
# Check if stop due to invalid bondlengths
valid_bondlengths = self.check_bondlengths()
if not valid_displace or not valid_z or not valid_bondlengths:
#self.r0 = self.atoms_prior.get_positions()
return True
def check_displacement(self):
valid_displace = True
if self.dmax_cov is not None:
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)
valid_displace = False
return valid_displace
def check_zcoord(self):
valid_z = True
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))
valid_z = False
return valid_z
def check_bondlengths(self):
valid_bondlengths = True
if self.blmin is not None or self.blmax is not None:
valid_bondlengths = check_valid_bondlengths(self.atoms, self.blmin, self.blmax)
if not valid_bondlengths:
#print('RELAXATION STOPPED on rank {} due to - invalid bondlengths.'.format(self.rk))
valid_bondlengths = False
return valid_bondlengths
candidate_operations/basic_mutations.py
View file @ 5c688cd1
......@@ -5,7 +5,7 @@ from ase.geometry import get_distances
from ase.visualize import view
from candidate_operations.candidate_generation import OffspringOperation
from candidate_operations.candidate_generation import CandidateGenerator
def pos_add_sphere(rattle_strength):
......@@ -28,7 +28,7 @@ def pos_add_sphere_shell(rmin, rmax):
np.cos(phi)])
return pos_add
class RattleMutation(OffspringOperation):
class RattleMutation(CandidateGenerator):
"""Class to perform rattle mutations on structures.
Parameters:
......@@ -52,13 +52,15 @@ class RattleMutation(OffspringOperation):
description: Name of the operation, which will be saved in
info-dict of structures, on which the operation is applied.
"""
def __init__(self, n_top, Nrattle=3, rattle_range=3, blmin=0.7, blmax=1.3,
description='RattleMutation'):
OffspringOperation.__init__(self, blmin=blmin, blmax=blmax)
def __init__(self, n_top, Nrattle=3, rattle_range=3, blmin=0.7, blmax=1.4,
force_all_bonds_valid=False, description='RattleMutation'):
CandidateGenerator.__init__(self, blmin=blmin, blmax=blmax,
force_all_bonds_valid=force_all_bonds_valid)
self.description = description
self.n_top = n_top
self.probability = Nrattle/n_top
self.rattle_range = rattle_range
self.force_all_bonds_valid = force_all_bonds_valid
def get_new_candidate(self, parents):
a = parents[0]
......@@ -83,23 +85,26 @@ class RattleMutation(OffspringOperation):
# Perform rattle operations in sequence.
for i in indices_to_rattle:
posi_0 = np.copy(a.positions[i])
for _ in range(100):
posi_0 = np.copy(a.positions[i])
# Perform rattle
pos_add = pos_add_sphere(self.rattle_range)
a.positions[i] += pos_add
# Check if rattle was valid
valid_bondlengths = self.check_valid_bondlengths(a)
valid_bondlengths = self.check_valid_bondlengths(a, indices=[i])
if not valid_bondlengths:
a.positions[i] = posi_0
else:
break
return a
class RattleMutation2(OffspringOperation):
if valid_bondlengths:
return a
else:
# If mutation is not successfull in supplied number
# of trials, return initial structure.
return atoms
class RattleMutation2(CandidateGenerator):
"""Class to perform rattle mutations on structures.
Parameters:
......@@ -119,12 +124,14 @@ class RattleMutation2(OffspringOperation):
description: Name of the operation, which will be saved in
info-dict of structures, on which the operation is applied.
"""
def __init__(self, n_top, Nrattle=3, blmin=0.7, blmax=1.3,
description='RattleMutation'):
OffspringOperation.__init__(self, blmin=blmin, blmax=blmax)
def __init__(self, n_top, Nrattle=3, blmin=0.7, blmax=1.4,
force_all_bonds_valid=False, description='RattleMutation'):
CandidateGenerator.__init__(self, blmin=blmin, blmax=blmax,
force_all_bonds_valid=force_all_bonds_valid)
self.description = description
self.n_top = n_top
self.probability = Nrattle/n_top
self.force_all_bonds_valid = force_all_bonds_valid
def get_new_candidate(self, parents):
"""Standardized candidate generation method for all mutation
......@@ -150,8 +157,8 @@ class RattleMutation2(OffspringOperation):
# Perform rattle operations in sequence.
for i in indices_to_rattle:
posi_0 = np.copy(a.positions[i])
for _ in range(100):
posi_0 = np.copy(a.positions[i])
j = np.random.randint(Nslab,Natoms)
# Perform rattle
......@@ -162,16 +169,21 @@ class RattleMutation2(OffspringOperation):
a.positions[i] = np.copy(a.positions[j]) + pos_add
# Check if rattle was valid
valid_bondlengths = self.check_valid_bondlengths(a)
valid_bondlengths = self.check_valid_bondlengths(a, indices=[i])
if not valid_bondlengths:
a.positions[i] = posi_0
else:
break
return a
if valid_bondlengths:
return a
else:
# If mutation is not successfull in supplied number
# of trials, return initial structure.
return atoms
class PermutationMutation(OffspringOperation):
class PermutationMutation(CandidateGenerator):
"""Class to perform permutation mutations on structures.
Parameters:
......@@ -192,12 +204,14 @@ class PermutationMutation(OffspringOperation):
info-dict of structures, on which the operation is applied.
"""
def __init__(self, n_top, Npermute=3, blmin=0.7, blmax=1.3,
description='PermutationMutation'):
OffspringOperation.__init__(self, blmin=blmin, blmax=blmax)
def __init__(self, n_top, Npermute=3, blmin=0.7, blmax=1.4,
force_all_bonds_valid=False, description='PermutationMutation'):
CandidateGenerator.__init__(self, blmin=blmin, blmax=blmax,
force_all_bonds_valid=force_all_bonds_valid)
self.description = description
self.n_top = n_top
self.probability = Npermute/n_top
self.force_all_bonds_valid = force_all_bonds_valid
def get_new_candidate(self, parents):
"""Standardized candidate generation method for all mutation
......@@ -239,12 +253,17 @@ class PermutationMutation(OffspringOperation):
a.positions[j_permute] = pos_i
# Check if rattle was valid
valid_bondlengths = self.check_valid_bondlengths(a)
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
else:
break
return a
if valid_bondlengths:
return a
else:
# If mutation is not successfull in supplied number
# of trials, return initial structure.
return atoms
candidate_operations/candidate_generation.py
View file @ 5c688cd1
......@@ -2,94 +2,78 @@ import numpy as np
from abc import ABC, abstractmethod
from ase.data import covalent_radii
from ase.geometry import get_distances
from ase import Atoms
from ase.visualize import view
class OffspringOperation(ABC):
"""Baseclass for mutation and crossover operations.
from utils import check_valid_bondlengths
class CandidateGenerator(ABC):
"""Baseclass for mutation and crossover operations as well
as the startgenerator.
Parameters:
blmin: The minimum allowed distance between atoms in units of
the covalent distance between atoms, where d_cov=r_cov_i+r_cov_j.
blmax: The maximum allowed distance, in units of the covalent
distance, from a single isolated atom to the closest atom. If
blmax=None, no constraint is enforced on isolated atoms.
force_all_bonds_valid: If True all bondlengths are forced to
be valid according to blmin/blmax. If False, only bondlengths
of atoms specified in bondlength checks during operations are
tested. The specified atoms are typically the ones changed
during operations. Default is False, as True might cause
problems with GOFEE, as GPR-relaxations and dual-steps might
result in structures that does not obey blmin/blmax.
"""
def __init__(self, blmin=0.7, blmax=1.3):
def __init__(self, blmin=0.7, blmax=1.4, force_all_bonds_valid=False):
self.blmin = blmin
self.blmax = blmax
self.force_all_bonds_valid = force_all_bonds_valid
self.description = 'Unspecified'
def check_valid_bondlengths(self, a, indices=None, check_too_close=True, check_isolated=True):
if self.force_all_bonds_valid:
# Check all bonds (mainly for testing)
return check_valid_bondlengths(a, self.blmin, self.blmax+0.1,
check_too_close=check_too_close,
check_isolated=check_isolated)
else:
# Check only specified ones
# (typically only for the atoms changed during operation)
return check_valid_bondlengths(a, self.blmin, self.blmax+0.1, indices=indices,
check_too_close=check_too_close,
check_isolated=check_isolated)
@abstractmethod
def get_new_candidate(self):
pass
def check_valid_bondlengths(self, a, indices=None):
bl = self.get_distances_as_fraction_of_covalent(a,indices)
# Filter away self interactions.
bl = bl[bl > 1e-6].reshape(bl.shape[0],bl.shape[1]-1)
# Check if atoms are too close
tc = np.any(bl < self.blmin)
# Check if there are isolated atoms
if self.blmax is not None:
isolated = np.any(np.all(bl > self.blmax, axis=1))
else:
isolated = False
is_valid = not tc and not isolated
return is_valid
def get_covalent_distance_from_atom_numbers(self, a, indices=None):
r_cov = np.array([covalent_radii[n] for n in a.get_atomic_numbers()])
if indices is None:
r_cov_sub = r_cov
else:
r_cov_sub = r_cov[indices]
cd_mat = r_cov_sub.reshape(-1,1) + r_cov.reshape(1,-1)
return cd_mat
def get_distances_as_fraction_of_covalent(self, a, indices=None, covalent_distances=None):
if indices is None:
indices = np.arange(len(a))
if covalent_distances is None:
cd = self.get_covalent_distance_from_atom_numbers(a, indices=indices)
def finalize(self, a, a0=None, successfull=True):
if successfull:
description = self.description
else:
cd = covalent_distances[indices,:]
_, d = get_distances(a[indices].positions,
a.positions,
cell=a.get_cell(),
pbc=a.get_pbc())
bl = d/cd
return bl
def finalize(self, a):
a.info['key_value_pairs'] = {'origin': self.description}
valid_bondlengths = self.check_valid_bondlengths(a)
assert valid_bondlengths, 'bondlengths are not valid'
description = 'failed ' + self.description
try:
a.info['key_value_pairs']['origin'] = description
except:
a.info['key_value_pairs'] = {'origin': description}
if self.force_all_bonds_valid:
# Check all bonds
valid_bondlengths = self.check_valid_bondlengths(a)
assert valid_bondlengths, 'bondlengths are not valid'
return a
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_sphere_shell(rmin, rmax):
# Rattle within a sphere
r = np.random.uniform(rmin**3, rmax**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)])
return pos_add
class CandidateGenerator():
"""Class to generate new candidates by applying candidate generation
operations drawn randomly from a list of possible operations
class OperationSelector():
"""Class to produce new candidates by applying one of the
candidate generation operations which is supplied in the
'operations'-list. The operations are drawn randomly according
to the 'probabilities'-list.
operations : "list" or "list of lists" of mutations/crossovers.
......@@ -128,33 +112,159 @@ class CandidateGenerator():
parents[0] = anew
return anew
def random_pos(box):
""" Returns a random position within the box
described by the input box. """
p0 = box[0].astype(float)
vspan = box[1]
r = np.random.random((1, len(vspan)))
pos = p0.copy()
for i in range(len(vspan)):
pos += vspan[i] * r[0, i]
return pos
class StartGenerator(CandidateGenerator):
""" Class used to generate random initial candidates.
The candidates are generated by iteratively adding in
one atom at a time within the box described.
Parameters:
slab: The atoms object describing the super cell to
optimize within. Can be an empty cell or a cell
containing the atoms of a slab.
stoichiometry: A list of atomic numbers for the atoms
that are placed on top of the slab (if one is present).
box_to_place_in: The box within which atoms are placed. The box
should be on the form [p0, vspan] where 'p0' is the position of
the box corner and 'vspan' is a matrix containing the three
spanning vectors.
blmin: The minimum allowed distance between atoms in units of
the covalent distance between atoms, where d_cov=r_cov_i+r_cov_j.
blmax: The maximum allowed distance, in units of the covalent
distance, from a single isolated atom to the closest atom. If
blmax=None, no constraint is enforced on isolated atoms.
cluster: If True atoms are required to be placed within
blmin*d_cov of one of the other atoms to be placed. If
False the atoms in the slab are also included.
"""
def __init__(self, slab, stoichiometry, box_to_place_in,
blmin=0.7, blmax=1.4, cluster=False, description='StartGenerator'):
CandidateGenerator.__init__(self, blmin=blmin, blmax=blmax)
self.slab = slab
self.stoichiometry = stoichiometry
self.box = box_to_place_in
self.cluster = cluster
self.description = description
def get_new_candidate(self, parents=None):
a = self.make_structure()
a = self.finalize(a)
return a
def make_structure(self):
""" Generates a new random structure """
Nslab = len(self.slab)
Ntop = len(self.stoichiometry)
num = np.random.permutation(self.stoichiometry)
for i_trials in range(1000):
a = self.slab.copy()
for i in range(Ntop):
pos_found = False
for _ in range(300):
# Place new atom
posi = random_pos(self.box)
a += Atoms([num[i]], posi.reshape(1,3))
# Check if position of new atom is valid
not_too_close = self.check_valid_bondlengths(a, indices=[Nslab+i],
check_too_close=True,
check_isolated=False)
if len(a) == 1: # The first atom
not_isolated = True
else:
if self.cluster: # Check isolation excluding slab atoms.
not_isolated = self.check_valid_bondlengths(a[Nslab:], indices=[Nslab+i],
check_too_close=False,
check_isolated=True)
else: # All atoms.
not_isolated = self.check_valid_bondlengths(a, indices=[Nslab+i],
check_too_close=False,
check_isolated=True)
valid_bondlengths = not_too_close and not_isolated
if not valid_bondlengths:
del a[-1]
else:
pos_found = True
break
if not pos_found:
break
if pos_found:
break
if i_trials == 999 and not pos_found:
raise RuntimeError('StartGenerator: No valid structure was produced in 1000 trials.')
else:
return a
if __name__ == '__main__':
from ase.io import read
from ase.visualize import view
from candidate_operations.basic_mutations import RattleMutation, RattleMutation2, PermutationMutation
print(0.7*2*covalent_radii[1], 1.3*2*covalent_radii[1])
np.random.seed(1)
np.random.seed(7)
#a = read('/home/mkb/DFT/gpLEA/Si3x3/ref/gm_unrelaxed_done.traj', index='0')
#a = read('si3x3.traj', index='0')
#a = read('c6h6.traj', index='0')
a = read('sn2o3.traj', index='0')
traj = read('c6h6_init.traj', index=':')
#a = read('sn2o3.traj', index='0')
#slab = read('slab_sn2o3.traj', index='0')
"""
stoichiometry = 6*[50] + 10*[8]
c = slab.get_cell()
c[2,2] = 3.3
p0 = np.array([0,0,14])
box = [p0, c]
rattle = RattleMutation(n_top=16, Nrattle=3)
sg = StartGenerator(slab, stoichiometry, box)
"""
a = traj[0]
rattle = RattleMutation(n_top=len(a), Nrattle=3, rattle_range=2)
rattle2 = RattleMutation2(n_top=16, Nrattle=0.1)
permut = PermutationMutation(n_top=16, Npermute=2)
candidategenerator = CandidateGenerator([0., 1., 0.], [rattle, rattle2, permut])
candidategenerator = OperationSelector([1], [rattle])