import itertools


def four_body_sum(spins):
    """Calculate four body term in the periodic lattice
    Input: 
    spins: spins configuration matrix.
    Output:
    sum of four body terms.
    """
    size = len(spins)
    Sum = 0
    for i,j in itertools.product(range(size), repeat=2):
        Sum += spins[i,j] * spins[i,(j+1)%size] * spins[(i+1)%size,j] \
                    * spins[(i+1)%size, (j+1)%size] # only consider the right and top part for each site
    return Sum