Different choices of initialization of the search should now work....

Different choices of initialization of the search should now work. Particularly structures/startgenerator and candidate_operations

candidate_operations/candidate_generation.py

@@ -89,7 +89,7 @@ class OperationSelector():
         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
+        assert element_count_operations == element_count_probabilities, 'the two lists must have the same shape'
         self.operations = operations
         self.rho = [np.cumsum(prob_list) for prob_list in probabilities]
 

gofee.py

@@ -8,11 +8,17 @@ from population import population
 from ase import Atoms
 from ase.io import read, write, Trajectory
 from ase.calculators.singlepoint import SinglePointCalculator
+from ase.calculators.calculator import FileIOCalculator
 
 from bfgslinesearch_zlim import BFGSLineSearch_zlim
 from bfgslinesearch_constrained import BFGSLineSearch_constrained
 from ase.ga.relax_attaches import VariansBreak
 from parallel_utils import split, parallel_function_eval
+
+from candidate_operations.candidate_generation import OperationSelector
+from candidate_operations.basic_mutations import RattleMutation
+
+
 from mpi4py import MPI
 world = MPI.COMM_WORLD
 
@@ -43,7 +49,23 @@ def relax(structure, calc, Fmax=0.05, steps_max=200, dmax_cov=None):
 class GOFEE():
     """GOFEE global structure search method.
         
-    structures: 
+    structures: Atoms-object, list of Atoms-objects or None.
+    In initial structures from which to start the sesarch.
+    If None, the startgenerator must be supplied.
+    If less than Ninit structures is supplied, the remaining
+    ones are generated using the startgenerator or by rattling
+    the supplied structures, depending on wether the
+    startgenerator is supplied.
+
+    gpr: The Gaussian Process Regression model used as the
+    surrogate model for the Potential energy surface.
+    
+    startgenerator: 
+
+    trajectory:
+
+    kappa: 
+    
     """
     def __init__(self, structures=None,
                  calc=None,
@@ -52,7 +74,7 @@ class GOFEE():
                  candidate_generator=None,
                  trajectory=None,
                  kappa=2,
-                 Neval=200,
+                 Niter=200,
                  Ninit=10,
                  dmax_cov=3.5,
                  Ncandidates=30,
@@ -93,10 +115,23 @@ class GOFEE():
         else:
             self.startgenerator = startgenerator
 
+        if startgenerator is not None:
+            self.n_to_optimize = len(self.startgenerator.stoichiometry)
+        else:
+            self.n_to_optimize = len(self.structures[0])
+            for constraint in self.structures[0].constraints:
+                if isinstance(constraint, FixAtoms):
+                    indices_fixed = constraint.get_indices()
+                    self.n_to_optimize -= len(indices_fixed)
+                    break
+            
         if candidate_generator is not None:
             self.candidate_generator = candidate_generator
         else:
-            assert False
+            rattle = RattleMutation(self.n_to_optimize,
+                                    Nrattle=3,
+                                    rattle_range=4)
+            self.candidate_generator = OperationSelector([1.0],[rattle])
 
         # Initialize population
         self.population = population(population_size=population_size, gpr=self.gpr, similarity2equal=0.9999)
@@ -106,7 +141,7 @@ class GOFEE():
         self.master = self.comm.rank == 0
 
         self.kappa = kappa
-        self.Neval = Neval
+        self.Niter = Niter
         self.Ninit = Ninit
         self.dmax_cov = dmax_cov
         self.Ncandidates = Ncandidates
@@ -133,30 +168,60 @@ class GOFEE():
         self.steps = 0
 
     def get_initial_structures(self):
-        # Initial structures
-
+        """Method to prepare the initial structures for the search.
+        
+        The method makes sure that there are atleast self.Ninit
+        initial structures.
+        These structures are first of all the potentially supplied
+        structures. If more structures are required, these are
+        generated using self.startgenerator (if supplied), otherwise
+        they are generated by heavily rattling the supplied structures.
+        """
+        # Collect potentially supplied structures and evaluate
+        # energies and forces if not present.
+        structures_init = []
         if self.structures is not None:
             for a in self.structures:
-                Ei = a.get_potential_energy()
-                Fi = a.get_forces()
-                self.gpr.memory.save_data([a])
-                self.trajectory.write(a)
-                self.population.add(a)
-        elif self.startgenerator is not None:
-            for i in range(int(np.ceil(self.Ninit/2))):
-                a = self.startgenerator.get_new_candidate()
+                calc = a.get_calculator()
+                if isinstance(calc, SinglePointCalculator):
+                    if 'energy' in calc.results and 'forces' in calc.results:
+                        # Add without evaluating.
+                        structures_init.append(a)
+                        self.write(a)
+                        continue
                 a = self.evaluate(a)
-                a2 = self.get_dualpoint(a) 
-                a2 = self.evaluate(a2)
-                self.gpr.memory.save_data([a,a2])
-                self.population.add(a2)
+                structures_init.append(a)
+
+        Nremaining = self.Ninit - len(structures_init)
+    
+        if Nremaining > 0 and self.startgenerator is None:
+            # Initialize rattle-mutation for all atoms.
+            rattle = RattleMutation(self.n_to_optimize,
+                                    Nrattle=self.n_to_optimize,
+                                    rattle_range=2)
+
+        # Generation of remaining initial-structures (up to self.Ninit).
+        for i in range(Nremaining):
+            if self.startgenerator is not None:
+                a = self.startgenerator.get_new_candidate()
             else:
-            assert False
+                # Perform two times rattle of all atoms.
+                a0 = structures_init[i % len(structures_init)]
+                a = rattle.get_new_candidate([a])
+                a = rattle.get_new_candidate([a])
+            a = self.evaluate(a)
+            structures_init.append(a)
+
+        # Potentially do a few relaxation steps.
+        ### missing code ###
+            
+        self.gpr.memory.save_data(structures_init)
+        self.population.add(structures_init)
                 
     def run(self):
         self.get_initial_structures()
 
-        while self.steps < self.Neval:
+        while self.steps < self.Niter:
             self.print_master('steps:', self.steps)
             self.train_surrogate()
             self.update_population()
@@ -315,10 +380,19 @@ class GOFEE():
     def evaluate(self, a):
         a = self.comm.bcast(a, root=0)
         
+        if isinstance(self.calc, FileIOCalculator):
+            if self.master:
+                a.set_calculator(self.calc)
+                E = a.get_potential_energy()
+                F = a.get_forces()
+                results = {'energy': E, 'forces': F}
+                calc_sp = SinglePointCalculator(a, **results)
+                a.set_calculator(calc_sp)
+            a = self.comm.bcast(a, root=0)
+        else:
             a.set_calculator(self.calc)
             E = a.get_potential_energy()
             F = a.get_forces()
-        
             results = {'energy': E, 'forces': F}
             calc_sp = SinglePointCalculator(a, **results)
             a.set_calculator(calc_sp)