update

candidate_operations/basic_mutations.py

+import numpy as np
+from abc import ABC, abstractmethod
+from ase.data import covalent_radii
+from ase.geometry import get_distances
+
+from ase.visualize import view
+
+from candidate_operations.candidate_generation import OffspringOperation
+
+
+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 RattleMutation(OffspringOperation):
+    """Class to perform rattle mutations on structures.
+    
+    Parameters:
+
+    n_top: The number of atoms to optimize. Specifically the
+    atoms with indices [-n_top:] are optimized.
+    
+    Nrattle: The average number of atoms to rattle.
+
+    rattle_range: The maximum distance within witch to rattle the
+    atoms. Atoms are rattled uniformly within a sphere of this
+    radius.
+
+    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.
+
+    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)
+        self.description = description
+        self.n_top = n_top
+        self.probability = Nrattle/n_top
+        self.rattle_range = rattle_range
+
+    def get_new_candidate(self, parents):
+        a = parents[0]
+        a = self.rattle(a)
+        a = self.finalize(a)
+        return a
+
+    def rattle(self, atoms):
+        """Standardized candidate generation method for all mutation
+        and crossover operations.
+        """
+        a = atoms.copy()
+        Natoms = len(a)
+        Nslab = Natoms - self.n_top
+
+        # Randomly select indices of atoms to permute - in random order.
+        indices_to_rattle = np.arange(Nslab,Natoms)[np.random.rand(self.n_top)
+                                                     < self.probability]
+        indices_to_rattle = np.random.permutation(indices_to_rattle)
+        if len(indices_to_rattle) == 0:
+            indices_to_rattle = [np.random.randint(Nslab,Natoms)]
+
+        # Perform rattle operations in sequence.
+        for i in indices_to_rattle:
+            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)
+                
+                if not valid_bondlengths:
+                    a.positions[i] = posi_0
+                else:
+                    break
+        return a
+
+class RattleMutation2(OffspringOperation):
+    """Class to perform rattle mutations on structures.
+    
+    Parameters:
+
+    n_top: The number of atoms to optimize. Specifically the
+    atoms with indices [-n_top:] are optimized.
+    
+    Nrattle: The average number of atoms to rattle.
+
+    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.
+
+    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)
+        self.description = description
+        self.n_top = n_top
+        self.probability = Nrattle/n_top
+
+    def get_new_candidate(self, parents):
+        """Standardized candidate generation method for all mutation
+        and crossover operations.
+        """
+        a = parents[0]
+        a = self.rattle(a)
+        a = self.finalize(a)
+        return a
+
+    def rattle(self, atoms):
+        a = atoms.copy()
+        Natoms = len(a)
+        Nslab = Natoms - self.n_top
+        num = a.numbers
+
+        # Randomly select indices of atoms to permute - in random order.
+        indices_to_rattle = np.arange(Nslab,Natoms)[np.random.rand(self.n_top)
+                                                     < self.probability]
+        indices_to_rattle = np.random.permutation(indices_to_rattle)
+        if len(indices_to_rattle) == 0:
+            indices_to_rattle = [np.random.randint(Nslab,Natoms)]
+
+        # Perform rattle operations in sequence.
+        for i in indices_to_rattle:
+            for _ in range(100):
+                posi_0 = np.copy(a.positions[i])
+                j = np.random.randint(Nslab,Natoms)
+
+                # Perform rattle
+                covalent_dist_ij = covalent_radii[num[i]] + covalent_radii[num[j]]
+                rmin = self.blmin * covalent_dist_ij
+                rmax = self.blmax * covalent_dist_ij
+                pos_add = pos_add_sphere_shell(rmin, rmax)
+                a.positions[i] = np.copy(a.positions[j]) + pos_add
+
+                # Check if rattle was valid
+                valid_bondlengths = self.check_valid_bondlengths(a)
+                
+                if not valid_bondlengths:
+                    a.positions[i] = posi_0
+                else:
+                    break
+        return a
+
+
+class PermutationMutation(OffspringOperation):
+    """Class to perform permutation mutations on structures.
+    
+    Parameters:
+
+    n_top: The number of atoms to optimize. Specifically the
+    atoms with indices [-n_top:] are optimized.
+    
+    Npermute: The average number of permutations to perform.
+    
+    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.
+    
+    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, Npermute=3, blmin=0.7, blmax=1.3,
+                 description='PermutationMutation'):
+        OffspringOperation.__init__(self, blmin=blmin, blmax=blmax)
+        self.description = description
+        self.n_top = n_top
+        self.probability = Npermute/n_top
+
+    def get_new_candidate(self, parents):
+        """Standardized candidate generation method for all mutation
+        and crossover operations.
+        """
+        a = parents[0]
+        a = self.mutate(a)
+        a = self.finalize(a)
+        return a
+
+    def mutate(self, atoms):
+        a = atoms.copy()
+        Natoms = len(a)
+        Nslab = Natoms - self.n_top
+        num = a.numbers
+
+        # Check if permutation mutation is applicable to structure.
+        num_unique_top = list(set(num[-self.n_top:]))
+        assert len(num_unique_top) > 1, 'Permutations with one atomic type is not valid'
+
+        # Randomly select indices of atoms to permute - in random order.
+        indices_to_permute = np.arange(Nslab,Natoms)[np.random.rand(self.n_top)
+                                                   < self.probability]
+        indices_to_permute = np.random.permutation(indices_to_permute)
+        if len(indices_to_permute) == 0:
+            indices_to_permute = [np.random.randint(Nslab,Natoms)]
+
+        # Perform permutations in sequence.
+        for i_permute 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)
+
+                # 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
+
+                # Check if rattle was valid
+                valid_bondlengths = self.check_valid_bondlengths(a)
+                
+                if not valid_bondlengths:
+                    a.positions[i_permute] = pos_i
+                    a.positions[j_permute] = pos_j
+                else:
+                    break
+        return a
+

candidate_operations/candidateGenerator.py

-import numpy as np
-
-
-
-class candidateGenerator():
-    def __init__(self):
-        pass
-
-
-
-
-class candidateOperation():
-    def __init__(self):
-        pass

candidate_operations/candidate_generation.py

+import numpy as np
+from abc import ABC, abstractmethod
+from ase.data import covalent_radii
+from ase.geometry import get_distances
+
+from ase.visualize import view
+
+class OffspringOperation(ABC):
+    """Baseclass for mutation and crossover operations.
+    """
+    def __init__(self, blmin=0.7, blmax=1.3):
+        self.blmin = blmin
+        self.blmax = blmax
+        self.description = 'Unspecified'
+
+    @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)
+        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'
+        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
+    
+    operations : "list" or "list of lists" of mutations/crossovers.
+
+    probabilities : probability for each of the mutations/crossovers
+    in operations. Must have the same dimensions as operations.
+    """
+    def __init__(self, probabilities, operations):
+        cond1 = isinstance(operations[0], list)
+        cond2 = isinstance(probabilities[0], list)
+        if not cond1 and not cond2:
+            operations = [operations]
+            probabilities = [probabilities]
+        element_count_operations = [len(op_list) for op_list in operations]
+        element_count_probabilities = [len(prob_list)
+                                       for prob_list in probabilities]
+        assert element_count_operations == element_count_probabilities
+        self.operations = operations
+        self.rho = [np.cumsum(prob_list) for prob_list in probabilities]
+
+    def __get_index__(self, rho):
+        """Draw from the cumulative probalility distribution, rho,
+        to return the index of which operation to use"""
+        v = np.random.random() * rho[-1]
+        for i in range(len(rho)):
+            if rho[i] > v:
+                return i
+        
+    def get_new_candidate(self, parents):
+        """Generate new candidate by applying a randomly drawn
+        operation on the structures. This is done successively for
+        each list of operations, if multiple are present.
+        """
+        for op_list, rho_list in zip(self.operations, self.rho):
+            to_use = self.__get_index__(rho_list)
+            anew = op_list[to_use].get_new_candidate(parents)
+            parents[0] = anew
+        return anew
+
+if __name__ == '__main__':
+    from ase.io import read
+    from ase.visualize import view
+
+    from candidate_operations.basic_mutations import RattleMutation, RattleMutation2, PermutationMutation
+    
+    np.random.seed(1)
+    
+    #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')
+    
+    rattle = RattleMutation(n_top=16, Nrattle=3)
+    rattle2 = RattleMutation2(n_top=16, Nrattle=0.1)
+    permut = PermutationMutation(n_top=16, Npermute=2)
+
+    candidategenerator = CandidateGenerator([0., 1., 0.], [rattle, rattle2, permut])
+    #candidategenerator = CandidateGenerator([[1],[1]], [[rattle2], [permut]])
+
+    traj = []
+    for i in range(100):
+        a0 = a.copy()
+        anew = candidategenerator.get_new_candidate([a0,a0])
+        traj += [a0, anew]
+
+    view(traj)
+    """
+    a_mut = rattle.get_new_candidate([a])
+    view([a,a_mut])
+    """

candidate_operations/candidate_generation.py~

+import numpy as np
+from abc import ABC, abstractmethod
+from ase.data import covalent_radii
+from ase.geometry import get_distances
+
+class OffspringCreator(ABC):
+    def __init__(self, blmin=0.7, blmax=1.3):
+        self.blmin = blmin
+        self.blmax = blmax
+        self.description = 'Unspecified'
+
+    @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)
+        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
+        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
+    
+class RattleMutation(OffspringCreator):
+    def __init__(self, n_top, Nrattle=3, rattle_range=3, blmin=0.7, blmax=1.3):
+        OffspringCreator.__init__(self, blmin=blmin, blmax=blmax)
+        self.description = 'RattleMutation'
+        self.n_top = n_top
+        self.probability = Nrattle/n_top
+        self.rattle_range = rattle_range
+
+    def get_new_candidate(self, parents):
+        a = parents[0]
+        a = self.rattle(a)
+        a = self.finalize(a)
+        return a
+
+    def rattle(self, atoms):
+        a = atoms.copy()
+        Natoms = len(a)
+        Nslab = Natoms - self.n_top
+
+        for i in range(Nslab,Natoms):
+            if np.random.random() < self.probability:
+                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, indices=[i])
+
+                    if not valid_bondlengths:
+                        a.positions[i] = posi_0
+                    else:
+                        break
+        return a
+                    
+
+class candidateOperation():
+    def __init__(self):
+        pass
+
+
+if __name__ == '__main__':
+    from ase.io import read
+    from ase.visualize import view
+    
+    a = read('/home/mkb/DFT/gpLEA/Si3x3/ref/gm_unrelaxed_done.traj', index='0')
+    
+    rattle = RattleMutation(n_top=16, Nrattle=3)
+
+    a_mut = rattle.get_new_candidate([a])
+    view([a,a_mut])

candidate_operations/offspring_creator.py

-"""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