Generates the required Sz operators, given the system parameter setup
       and the operator dictionary

def spin_z_op(param, oper):
    """Generates the required Sz operators, given the system parameter setup
       and the operator dictionary"""
    slaves = param['slaves']
    oper['Sz'] = np.array([spin_z(slaves, spin) for spin in range(slaves)])
    oper['Sz+1/2'] = oper['Sz'] + 0.5*np.eye(2**slaves)
    oper['sumSz2'] = oper['Sz'].sum(axis=0)**2  # because Sz is diagonal
    Sz_mat_shape = oper['Sz'].reshape(param['orbitals'], 2, 2**slaves, 2**slaves)
    oper['sumSz-sp2'] = (Sz_mat_shape.sum(axis=1)**2).sum(axis=0)
    oper['sumSz-or2'] = (Sz_mat_shape.sum(axis=0)**2).sum(axis=0)