code: def symmetrically_remove_atoms(self, indices):
  
          
  
          slabcopy = SpacegroupAnalyzer(self.copy()).get_symmetrized_structure()
          points = [slabcopy[i].frac_coords for i in indices]
          removal_list = []
          
          for pt in points:
              # Get the index of the original site on top
              cart_point = slabcopy.lattice.get_cartesian_coords(pt)
              dist = [site.distance_from_point(cart_point) for site in slabcopy]
              site1 = dist.index(min(dist))
  
              # Get the index of the corresponding site at the bottom
              for i, eq_sites in enumerate(slabcopy.equivalent_sites):
                  if slabcopy[site1] in eq_sites:
                      eq_indices = slabcopy.equivalent_indices[i]
                      break
              i1 = eq_indices[eq_sites.index(slabcopy[site1])]
              
              for i2 in eq_indices:
                  if i2 == i1:
                      continue
                  if slabcopy[i2].frac_coords[2] == slabcopy[i1].frac_coords[2]:
                      continue
                  # Test site remove to see if it results in symmetric slab
                  s = self.copy()
                  s.remove_sites([i1, i2])
                  if s.is_symmetric():
                      removal_list.extend([i1, i2])
                      break
                      
          # If expected, 2 atoms are removed per index
          if len(removal_list) == 2*len(indices):
              self.remove_sites(removal_list)
          else:
              warnings.warn("Equivalent sites could not be found for removal for all indices. Surface unchanged.")
docs: Class method for removing sites corresponding to a list of indices.
              Will remove the corresponding site on the other side of the
              slab to maintain equivalent surfaces.
  
          Arg:
              indices ([indices]): The indices of the sites
                  in the slab to remove.