Using invariant theory, new invariant cubature formulas over a unit cube are given by imposing a group structure on the formulas. Cools and Haegemans [Computing 40, 139-146 (1988)] constructed invariant cubature formulas over a unit square. Since there exists a problem in directly extending their ideas over the unit square which were obtained by using a concept of good integrity basis to some constructions of invariant cubature formulas over the unit cube, a Reynold operator will be used to obtain new invariant cubature formulas over the unit cube. In order to practically find integration nodes and weights for the cubature formulas, it is required to solve a system of nonlinear equations. With an IMSL subroutine DUNLSF which is used for solutions of the system of nonlinear equations, we shall give integration nodes for the new invariant cubature formulas over the unit cube depending on each degree of polynomial precision.