An approach with statistical mechanics and a unified molecular theory of non-linear viscoelasticity with constraints of Gaussian chain entanglement for polymer melts were proposed. A multimode model structure for a single polymer chain with tail segments and N reversible entanglement sites on the test polymer chain was developed. The probability distribution function of the end-to-end vector for a single polymer chain at entangled state and the viscoelastic free energy of deformation for polymer melts were calculated. Four types of stress-strain relationship and memory function were derived from this theory. The above theoretical relationships were verified by experimental data of PS (polystyrene) and LDPE (low density polyethylene) melts.