structures are saved without the surface now, which makes some things easier;...

structures are saved without the surface now, which makes some things easier; side chain generation is separated from the reduction part; first modifications of the reduction part to work with arbitrary representations and selection algorithms

new/codered.py

@@ -4,6 +4,7 @@ from .adAtom import AdAtom
 from .sideChainGen import SideChain
 from .boundaryCondition import *
 from .reduction import *
+import .representations as representations
 from .utils import get_periodic_images,find_NNs
 from .output import printProgressBar
 from .errors import *
@@ -140,88 +141,19 @@ class CodeRed:
 
     def get_structures(self, reduced=False):
         """
-        This method returns a copy of the current structures array.
+        This method returns a copy of the current structures array including the surface.
         The surface can be replaced by the reduced surface by setting reduced=True.
         """
         structures = []
         for struct in self.structures:
-            copy = struct.copy()
             if reduced:
-                structures.append(self.surfaceRed.copy()+copy[len(self.surface):])
+                structures.append(self.surfaceRed.copy()+struct.copy())
             else:
-                structures.append(copy)
+                structures.append(self.surface.copy()+struct.copy())
 
         return structures
 
 
-    def make_sidechains(self):
-        """
-        This function calculates all side chain positions and adds the side chains to the 
-        current structures. For all side chains agreement with the minimum distance condition
-        within the side chain is checked, and if the condition is not met for any side chain
-        the corresponding structure is removed from the list of structures.
-        """
-        # get the corresponding indices of the connect_to attribute of all side chains
-        # (so it doesn't need to be done again and again in the loop)
-        sc_connect_indices = []
-        for sc in self.sideChains:
-            if isinstance(sc.connect_to, AdAtom):
-                if sc.connect_to in self.adAtoms:
-                    sc_connect_indices.append(self.adAtoms.index(sc.connect_to))
-                else:
-                    print('ERROR! Unable to find the AdAtom index connected with this side chain!')
-                    
-            elif isinstance(sc.connect_to, int):
-                if sc.connect_to < self.current_index:
-                    # convert to positive int, if negative
-                    sc_connect_indices.append(sc.connect_to % self.current_index)
-                else:
-                    print('ERROR! Unable to find the AdAtom index connected with this side chain!')
-
-            elif isinstance(sc.connect_to, str):
-                index = -1
-                for i in range(len(self.adAtoms)):
-                    if self.adAtoms[i].label==sc.connect_to:
-                        index = i
-                        break
-                    
-                if index!=-1:
-                    sc_connect_indices.append(index)
-                else:
-                    print('ERROR! Unable to find the AdAtom index connected with this side chain!')
-
-            else:
-                print('ERROR! Unable to recognise connect_to identifier for this side chain!')
-
-        # get current structures and add side chains
-        structures = self.get_structures()
-        new_structures = []
-
-        # create progress bar
-        progress = 0
-        printProgressBar(progress, len(structures), prefix=' Progress: ')
-        for struct in structures:
-            # calculate and add side chains
-            success = True
-            for i in range(len(self.sideChains)):
-                try:
-                    sc = self.sideChains[i].calc_positions(len(self.surface)+sc_connect_indices[i], struct)
-                    struct = struct + sc
-                except DistanceError:
-                    success = False
-                    break
-
-            # add to structure list
-            if success:
-                new_structures.append(struct)
-
-            # update progress bar
-            progress += 1
-            printProgressBar(progress, len(structures), prefix=' Progress: ')
-
-        return new_structures
-
-
     def code(self):
         """
         This is the core method to calculate the structures that comply
@@ -335,7 +267,7 @@ class CodeRed:
         for comb in self.combinations:
             # create atoms object from surface and current combo
             ad_structure = self.surface.copy() + Atoms(species, comb)
-            structures.append(ad_structure)
+            structures.append(ad_structure[-self.current_index:])
 
         # save list of structures
         self.structures = structures
@@ -343,47 +275,61 @@ class CodeRed:
         return
 
 
-    def red(self, red_type=None, reduced_surface=True, **reduction_params):
+    def red(self, red_type=None, surface='reduced', representation="positions", representation_params={}, selection="symmetry", selection_params={}, **reduction_params):
         #TODO: process reduction_params to set up and execute the right reduction technique
         #      e.g. symmetry reduction, mesh projection reduction, soap distance reduction, ...
 
         # get structures for reduction process
-        if reduced_surface:
+        if surface=='reduced' or surface==False:
             structures = self.get_structures(reduced=True)
-        else:
+        elif surface=='full' or surface==True:
             structures = self.get_structures()
-
-        # perform reduction according to the set reduction type
-        if red_type=='symmetry':
-            # generate informative output
-            print("Applying symmetry reduction...")
-            
-            # perform symmetry reduction
-            structures_red, red_array = symmetryReduction(structures, **reduction_params)
-            structures_new = []
-            for i in range(len(self.structures)):
-                if red_array[i]:
-                    structures_new.append(self.structures[i])
-
-        elif red_type=='grid':
-            # generate informative output
-            print("Applying grid based coarse graining...")
-
-            # perform grid based reduction
-            structures_red, structures_cg, red_array = gridCoarseGraining(structures, **reduction_params)
-            structures_new = []
-            for i in range(len(self.structures)):
-                if red_array[i]:
-                    structures_new.append(self.structures[i])
-
-        elif red_type=='sidechains':
-            # generate informative output
-            print("Generating side chains and corresponding minimum distance reduction...")
-
-            # generate structures including side chains
-            structures_new = self.make_sidechains()
-
-        # save new structures
+        elif surface=='none' or surface==None:
+            structures = self.structures.copy()
+
+        # transform structures according to representation setting
+        if representation in ["positions", "pos", "coordinates"]:
+            structures_repr = representations.positions(structures, **representation_params)
+        elif representation in ["box", "boxes", "boxgrid", "grid"]:
+            structures_repr = representations.boxgrid(structures, **representation_params)
+        else:
+            try:
+                structures_repr = representation(structures, **representation_params)
+            except TypeError:
+                print("This is not a valid representation!")
+                raise
+
+        ## perform reduction according to the set reduction type
+        #if red_type=='symmetry':
+        #    # generate informative output
+        #    print("Applying symmetry reduction...")
+        #    
+        #    # perform symmetry reduction
+        #    structures_red, red_array = symmetryReduction(structures, **reduction_params)
+        #    structures_new = []
+        #    for i in range(len(self.structures)):
+        #        if red_array[i]:
+        #            structures_new.append(self.structures[i])
+        #
+        #elif red_type=='grid':
+        #    # generate informative output
+        #    print("Applying grid based coarse graining...")
+        #
+        #    # perform grid based reduction
+        #    structures_red, structures_cg, red_array = gridCoarseGraining(structures, **reduction_params)
+        #    structures_new = []
+        #    for i in range(len(self.structures)):
+        #        if red_array[i]:
+        #            structures_new.append(self.structures[i])
+        #
+        #elif red_type=='sidechains':
+        #    # generate informative output
+        #    print("Generating side chains and corresponding minimum distance reduction...")
+        #
+        #    # generate structures including side chains
+        #    structures_new = self.make_sidechains()
+
+        # save reduced structure list
         self.structures = structures_new
         
         # generate informative output
@@ -392,6 +338,104 @@ class CodeRed:
         return
 
 
-    def reset(self):
+    def make_sidechains(self):
+        """
+        This function calculates all side chain positions and adds the side chains to the 
+        current structures. For all side chains agreement with the minimum distance condition
+        within the side chain is checked, and if the condition is not met for any side chain
+        the corresponding structure is removed from the list of structures.
+        """
+        # get the corresponding indices of the connect_to attribute of all side chains
+        # (so it doesn't need to be done again and again in the loop)
+        sc_connect_indices = []
+        for sc in self.sideChains:
+            if isinstance(sc.connect_to, AdAtom):
+                if sc.connect_to in self.adAtoms:
+                    sc_connect_indices.append(self.adAtoms.index(sc.connect_to))
+                else:
+                    print('ERROR! Unable to find the AdAtom index connected with this side chain!')
+                    
+            elif isinstance(sc.connect_to, int):
+                if sc.connect_to < self.current_index:
+                    # convert to positive int, if negative
+                    sc_connect_indices.append(sc.connect_to % self.current_index)
+                else:
+                    print('ERROR! Unable to find the AdAtom index connected with this side chain!')
+
+            elif isinstance(sc.connect_to, str):
+                index = -1
+                for i in range(len(self.adAtoms)):
+                    if self.adAtoms[i].label==sc.connect_to:
+                        index = i
+                        break
+                    
+                if index!=-1:
+                    sc_connect_indices.append(index)
+                else:
+                    print('ERROR! Unable to find the AdAtom index connected with this side chain!')
+
+            else:
+                print('ERROR! Unable to recognise connect_to identifier for this side chain!')
+
+        # get current structures (including full surface) and add side chains
+        structures = self.get_structures()
+        new_structures = []
+
+        # informative output
+        print("Generating side chains and checking for minimum distances...")
+
+        # create progress bar
+        progress = 0
+        printProgressBar(progress, len(structures), prefix=' Progress: ')
+        for struct in structures:
+            # calculate and add side chains
+            success = True
+            for i in range(len(self.sideChains)):
+                try:
+                    sc = self.sideChains[i].calc_positions(len(self.surface)+sc_connect_indices[i], struct)
+                    struct = struct + sc
+                except DistanceError:
+                    success = False
+                    break
+
+            # add to structure list (without the surface)
+            if success:
+                new_structures.append(struct[len(self.surface):])
+
+            # update progress bar
+            progress += 1
+            printProgressBar(progress, len(structures), prefix=' Progress: ')
+
+        self.structures = new_structures
+        print("Done! Remaining structures: {:d}".format(len(self.structures)))
+
+
+    def reset_code(self):
         self.current_index = 0
         self.combinations = []
+
+
+    def reset_red(self):
+        # make structures from combinations
+        print("Recreating structures from combinations...", end="\r")
+
+        # assemble list of added species
+        species = ''
+        for i in range(self.current_index):
+            species += self.adAtoms[i].species
+
+        # create list of structures
+        structures = []
+        for comb in self.combinations:
+            # create atoms object from surface and current combo
+            ad_structure = self.surface.copy() + Atoms(species, comb)
+            structures.append(ad_structure[-self.current_index:])
+
+        # save list of structures
+        self.structures = structures
+        print("Recreating structures from combinations... Done!")
+
+
+    def reset(self):
+        self.reset_code()
+        self.structures = []

new/representations.py

+import numpy as np
+from .output import printProgressBar
+
+
+# TODO: add progress bars?
+
+def positions(structures, cartesian=False):
+    # get all positions in the input
+    representation = []
+    for struct in structures:
+        struct.wrap()
+        if cartesian:
+            representation.append(struct.get_positions())
+        else:
+            representation.append(struct.get_scaled_positions())
+    
+    return representation
+
+
+def boxgrid(structures, n_boxes_a=1, n_boxes_b=1, n_boxes_c=1, cartesian=False):
+    """
+    This representation splits the bounding box of the structures into a grid of
+    boxes and represents each position by the center of the sub-box it is in.
+    The number of sub-boxes in each direction can be set via the input parameters.
+    If cartesian is set to True, the boxgrid will be created along the cartesian
+    axes instead of the unit cell axes.
+    """
+    
+    # calculate the bounding box
+    positions = np.array(positions(structures, cartesian))
+    grid_amin = np.amin(positions[:,:,0])
+    grid_amax = np.amax(positions[:,:,0])
+    grid_bmin = np.amin(positions[:,:,1])
+    grid_bmax = np.amax(positions[:,:,1])
+    grid_cmin = np.amin(positions[:,:,2])
+    grid_cmax = np.amax(positions[:,:,2])
+
+    # create the sub-box center grid
+    ga, gwa = np.linspace(grid_amin, grid_amax, n_boxes_a, endpoint=False, retstep=True)
+    gb, gwb = np.linspace(grid_bmin, grid_bmax, n_boxes_b, endpoint=False, retstep=True)
+    gc, gwc = np.linspace(grid_cmin, grid_cmax, n_boxes_c, endpoint=False, retstep=True)
+
+    grid_aa, grid_bb, grid_cc = np.meshgrid(ga+0.5*gwa, gb+0.5*gwb, gc+0.5*gwc, indexing='ij')
+    grid = np.array([grid_aa, grid_bb, grid_cc]).reshape(3,len(ga)*len(gb)*len(gc)).T
+
+    # project positions on the grid
+    representation = np.copy(positions).reshape(-1,3)
+    for i in range(len(representation)):
+        representation[i] = grid[np.argmin(np.linalg.norm(grid-representation[i],axis=-1))]
+    
+    return representation.reshape(positions.shape)