def get_pdm_terms(site_index, n, adj_sites=4, shift=0):
    """
    inputs:
        site_index (int): center site index
        n (int): number of sites in lattice (for enforcing periodic boundary conditions)
        adj_site (int): how many adjacent sites to collect pop and coherence values from
        shift (int): return row indices shifted by this value (when collecting from matrix w/ compound index
    returns:
        ind_list (list of ints): all site indices to collect
        shift_ind_list (list of ints): site indices shifted down i*n rows to collect the correct sample
        coh_list (list of ints): coherences in fragment with central site
        shift_coh_list (list of ints): coherences shifted down i*n rows to select particular sample
    """

    # if a term in sites is out of bounds, subtract bound length
    ind_list, shift_ind_list = [site_index], [site_index+shift]
    coh_list, shift_coh_list = [], []
    #print("Site: ", site_index)
    for ind in range(site_index - adj_sites, site_index + adj_sites + 1):
        if ind != site_index: # we've already add the target site population to ind_list and shift_ind_list
            if ind < 0: ind += n
            elif ind >= n: ind -= n
            else: pass
            #print(ind)
            ind_list.append(ind)
            shift_ind_list.append(ind+shift) # shift down to specific matrix in set we are selecting
            coh_list.append(ind)
            shift_coh_list.append(ind+shift)

    return ind_list, shift_ind_list, coh_list, shift_coh_list