Methods of implementation of classical molecular dynamics simulations of moderate size molecule vibrational energy relaxation and analysis of their results are proposed. Two different approaches are considered. The first is concerned with modeling a real nonequilibrium cooling process for the excited molecule in a solvent initially at equilibrium. In addition to the solute total, kinetic, and potential energy evolution, that define the character of the process and the rate constant or relaxation time, a great deal of important information is provided by a normal mode specific analysis of the process. Expressions for the decay of the normal mode energies, the work done by particular modes, and the vibration-rotation interaction are presented. The second approach is based on a simulation of a solute-solvent system under equilibrium conditions. In the framework of linear nonequilibrium statistical thermodynamics and normal mode representation of the solute several expressions for the rate constant are derived. In initial form, they are represented by integrals of the time correlation functions of the capacities of the solute-solvent interaction atomic or normal mode forces and include the solute heat capacity. After some approximations, which are adequate for specific cases, these expressions are transformed to combinations of those for individual oscillators with force-force time correlation functions. As an attempt to consider a strongly nonequilibrium situation we consider a two-temperature model and discuss the reason why the rate constant can be independent on the solute energy or temperature. Expressions for investigation of the energy redistribution in the solvent are derived in two forms. One of them is given in the usual form of a heat transfer equation with the source term describing the energy flux from the excited solute. The other form describes the energy redistribution in the solvent in terms of capacity time correlation functions and can be more convenient if memory effects and spatial dispersion play an important role in energy redistribution in the solvent.