def state_hamiltonian(state):
    """
    Given a spin state, sum through each nearest neighbor pair to find the
    energy of the spin hamiltonian.
    
    A factor of 1/2 is used to adjust for double counting
    """
    N = state.shape[0]
    energy = 0
    for ii in range(N):
        for jj in range(N):
            
            # Sum through sites
            site = state[ii, jj]
            
            #Sum nearest neighbor bonds, using modulo N to enforce periodic boundaries
            nn = state[(ii + 1) % N, jj] + state[(ii - 1) % N, jj] + \
                 state[ii, (jj + 1) % N] + state[ii, (jj - 1) % N]
                            
            # Update energy, eliminate double counting
            energy += -0.5*site*nn
    return energy