code: def spin_z_op(param, oper):
      
      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)
docs: Generates the required Sz operators, given the system parameter setup
         and the operator dictionary