added grid based coarse graining reduction method

new/codered.py

@@ -2,7 +2,7 @@ import numpy as np
 from ase import Atoms
 from .adAtom import AdAtom
 from .boundaryCondition import *
-from .reduction import symmetryReduction
+from .reduction import *
 
 
 class CodeRed:
@@ -187,35 +187,37 @@ class CodeRed:
         return
 
 
-    def red(self, **reduction_params):
+    def red(self, red_type=None, reduced_surface=True, **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' in reduction_params and reduction_params['reduced_surface']:
+        if reduced_surface:
             structures = self.get_structures(reduced=True)
         else:
             structures = self.get_structures()
 
         # perform reduction according to the set reduction type
-        if reduction_params['red_type']=='symmetry':
+        if red_type=='symmetry':
             # generate informative output
             print("Applying symmetry reduction...")
-            # extract possible extra parameters
-            symred_params = {}
-            try:
-                symred_params['tolerance'] = reduction_params['tolerance']
-            except KeyError:
-                pass
-            try:
-                symred_params['to_primitive'] = reduction_params['to_primitive']
-            except KeyError:
-                pass
             
             # perform symmetry reduction
-            structures_red, red_array = symmetryReduction(structures, **symred_params)
-            self.combinations = self.combinations[red_array]
-            
+            structures_red, red_array = symmetryReduction(structures, **reduction_params)
+
+        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)
+
+        else:
+            red_array = np.ones(len(structures), dtype='bool')
+
+        # reduce combinations accordingly
+        self.combinations = self.combinations[red_array]
+        
         # generate informative output
         print("Done! Reduced to {:d} structures in total.".format(len(self.combinations)))
 

new/reduction.py

@@ -4,6 +4,11 @@ from .output import printProgressBar
 
 
 def symmetryReduction(structures, tolerance=0.001, to_primitive=False):
+    """
+    This method checks which of the given structures have a symmetrically equivalent
+    structure in the list by using ase's symmetry equivalence check.
+    """
+    
     # set parameters for the symmetry check function
     secheck = SymmetryEquivalenceCheck(stol=tolerance, to_primitive=to_primitive)
 
@@ -14,14 +19,70 @@ def symmetryReduction(structures, tolerance=0.001, to_primitive=False):
     symequiv = np.zeros(n_struct, dtype='bool')
 
     for i in range(n_struct):
+        symequiv[i] = secheck.compare(structures[i],structures[i+1:])
         printProgressBar(i+1, n_struct, prefix=' Progress: ')
-        # only continue if there was not found an equivalent structure yet
-        if not symequiv[i]:
-            for j in range(i+1, n_struct):
-                # only continue if there was not found an equivalent structure yet
-                if not symequiv[j]:
-                    # check if structures i and j are equivalent
-                    symequiv[j] = secheck.compare(structures[i],structures[j])
 
     # return the structures without an equivalent structure and the decision array
     return np.array(structures)[np.logical_not(symequiv)], np.logical_not(symequiv)
+
+
+def gridCoarseGraining(structures, ngrains_x=1, ngrains_y=1, ngrains_z=1):
+    """
+    This method projects the positions of the given structures on a coarse mesh
+    and then checks for equivalence of structures on that coarse mesh, but ignoring
+    symmetry equivalence.
+    (All the structures have to have the same unit cell)
+    """
+
+    n_struct = len(structures)
+    equiv = np.zeros(n_struct, dtype='bool')
+    structures_cg = []
+
+    # calculate boundaries for coarse graining grid
+    grid_xmin = 0.5
+    grid_xmax = 0.5
+    grid_ymin = 0.5
+    grid_ymax = 0.5
+    grid_zmin = 0.5
+    grid_zmax = 0.5
+    
+    for i in range(n_struct):
+        positions = structures[i].get_scaled_positions()
+        grid_xmin = np.amin([grid_xmin, np.amin(positions[:,0])])
+        grid_xmax = np.amax([grid_xmax, np.amax(positions[:,0])])
+        grid_ymin = np.amin([grid_ymin, np.amin(positions[:,1])])
+        grid_ymax = np.amax([grid_ymax, np.amax(positions[:,1])])
+        grid_zmin = np.amin([grid_zmin, np.amin(positions[:,2])])
+        grid_zmax = np.amax([grid_zmax, np.amax(positions[:,2])])
+
+    # create coarse graining grid
+    gx, gwx = np.linspace(grid_xmin, grid_xmax, ngrains_x+1, endpoint=True, retstep=True)
+    gy, gwy = np.linspace(grid_ymin, grid_ymax, ngrains_y+1, endpoint=True, retstep=True)
+    gz, gwz = np.linspace(grid_zmin, grid_zmax, ngrains_z+1, endpoint=True, retstep=True)
+
+    grid_xx, grid_yy, grid_zz = np.meshgrid(gx, gy, gz, indexing='ij')
+    grid = np.array([grid_xx, grid_yy, grid_zz]).reshape(3,len(gx)*len(gy)*len(gz)).T
+    
+    for i in range(n_struct):
+        # copy structure to not modify the original one
+        struct_cg = structures[i].copy()
+
+        # transform structure, if necessary
+        struct_cg.wrap()
+        
+        # project positions on the grid
+        positions = struct_cg.get_scaled_positions()
+        positions_transformed = positions[:,np.newaxis,:].repeat(len(grid),axis=1)
+        grid_transformed = grid[np.newaxis,:,:].repeat(len(positions),axis=0)
+        grid_dist = np.linalg.norm(positions_transformed-grid_transformed, axis=-1) # distances to grid positions
+
+        struct_cg.set_scaled_positions(grid[grid_dist.argmin(axis=1)])
+        
+        # test if already included in coarse grained structures
+        equiv[i] = struct_cg in structures_cg
+        if not equiv[i]:
+            # add to the list of (unique) coarse grained structures
+            structures_cg.append(struct_cg)
+
+    # return non-equivalent structures (original and cg positions) and the decision array
+    return np.array(structures)[np.logical_not(equiv)], structures_cg, np.logical_not(equiv)