The simplest quantum approach that offers the anharmonic treatment of the molecular vibrations is the secondorder perturbation theory, usually denoted as VPT2. In this approach, harmonic, cubic (∝ V 3 ξ 3 ) and semidiagonal quartic force fields (∝ V 4 ξ 4 ) are taken into consideration. Since the harmonic potential is symmetric and the cubic part is antisymmetric with respect to the displacements from the equilibrium values, the perturbative first order (PT1) energy correction with the cubic force field is zero. However, the quartic semi-diagonal part of the field with the terms of the form V 4,ij ξ 2 i ξ 2 j , where i and j enumerate displacements along different modes, gives nonzero corrections at the PT1 level. The second order corrections (PT2) are thus being performed only for the cubic part of the potential. Such a procedure can be automated and is already implemented in different software packages for quantum chemistry. However, with the increasing number of atoms N in a molecule the computational cost of VPT2 increases dramatically due to the growing amount of cubic force field elements, which scales as N 3 . For large molecules MD simulations may become computationally more effective. Thus it is important to develop procedures for computation of molecular vibrational properties at low temperatures using classical MD simulations.