Extraction of a feature for all atoms within a 1×1×1 computational cell at the center of a 3×3×3 supercell

Question

Context:

To account for periodic boundary conditions, feature extraction is performed for all atoms within a 1×1×1 computational cell at the centre of a 3×3×3 supercell.

Actually, I want to extract Anion framework coordination (.py file attached here) for the computational unit cell at the super cell's centre (3x3x3).

I did it for all the anions that are in the Supercell (3x3x3). But I am unable to do it for only anions that are in the central unit cell of the Supercell.

In case you would like to download the original unit cell structure file (.cif), and the .py file (also pasted below in the problem sub-heading) here is a Google Drive Link

Problem

In the below code, the result for AFC I am getting is 7.445, as it is calculating for all the anions. But I want only anions that are in the central unit cell, and the answer will be 12 in this case. Here anions for which AFC is calculated should be part of a unit cell, but when calculating coordinate, that anion can be a part of a supercell.

from pymatgen.io.cif import CifParser
import numpy as np

def calculate_coordination(structure, anion_index, rij0, processed_indices):
    """
    Calculates the coordination of a given anion in a structure

    Args:
        structure (Structure): Pymatgen Structure object
        anion_index (int): Index of the anion atom
        rij0 (float): Minimum distance to consider for coordination calculation
        processed_indices (list): List of indices of processed anions

    Returns:
        int: coordination number of the anion with respect to other anions in the structure
    """
    coordination = 0
    anion = structure[anion_index]
    for i in range(len(structure)):
        if i == anion_index or structure[i].specie.X < anion.specie.X or i in processed_indices:
            continue
        if structure[i].specie.X == anion.specie.X:
            distance = np.linalg.norm(structure[i].coords - anion.coords)
            if rij0 <= distance <= rij0 + 1:
                coordination += 1
                processed_indices.append(i)
    return coordination
def calculate_anion_framework_coordination(filepath, supercell_dim):
    """
    Calculates the average coordination number of anions in the 1x1x1 supercell at the center of a 3x3x3 supercell

    Args:
        filepath (str): Path to cif file
        supercell_dim (list of 3 ints): Dimensions of the supercell in the format [x, y, z]

    Returns:
        float: Average coordination number of anions in the 1x1x1 supercell at the center of a 3x3x3 supercell
    """
    
    # Read cif file
    parser = CifParser(filepath)
    structure = parser.get_structures()[0]

    # Print the number of atoms in the original unit cell
    print("Number of atoms in original unit cell:", len(structure))

    # Create 3x3x3 supercell
    structure.make_supercell(supercell_dim)

    # Print the number of atoms in the supercell
    print("Number of atoms in supercell:", len(structure))

    # Find atom with highest electronegativity
    electronegativity = [atom.specie.X for atom in structure]
    anion_index = np.argmax(electronegativity)
    print(anion_index)

    # Get the center of the 3x3x3 supercell
    supercell_center = np.mean(structure.cart_coords, axis=0)

    # Get the indices of the anion atoms in the 3x3x3 supercell
    anion_indices = [i for i, site in enumerate(structure) if site.specie.X == structure[anion_index].specie.X]
    total_number_of_anions = len(anion_indices)
    print("total_number_of_anions:", total_number_of_anions)

    # Calculate the coordination of the anions
    total_coordination = 0
    num_anions = 0
    
    for i in anion_indices:
        anion = structure[i]
        print(anion)
        distances = []
        for j in anion_indices:
            if i == j:
                continue
            distance = np.linalg.norm(structure[j].coords - anion.coords)
            distances.append(distance)
    

        rij0 = min(distances)
        print(f"rij0 for anion {i+1}: {rij0}")
        rij1 = rij0 + 1
        print(f"Distances for which coordination is calculated for anion {i+1}: [{rij0:.3f}, {rij1:.3f}]")
        coordination = calculate_coordination(structure, i, rij0, [])
        print(coordination)
        total_coordination += coordination
        total_number_of_anions = len(anion_indices)
    afc = total_coordination / total_number_of_anions
    return afc
        
      
afc = calculate_anion_framework_coordination("/Users/rahulsharma/Desktop/Li3P.cif", [3,3,3])
print("Anion framework coordination:", afc)