The theoretical study has been performed to refine the procedure for calculations of Gibbs free energy with a relative accuracy of less than 1 kcal/mol. Three benchmark intermolecular complexes are examined via several quantum-chemical methods, including the second-order Moller-Plesset perturbation (MP2), coupled cluster (CCSD(T)), and density functional (BLYP, B3LYP) theories augmented by Dunnings correlation-consistent basis sets. The effects of electron correlation, basis set size, and anharmonicity are systematically analyzed, and the results are compared with available experimental data. The results of the calculations suggest that experimental accuracy can be reached only by extrapolation of MP2 and CCSD(T) total energies to the complete basis set. The contribution of anharmonicity to the zero point energy and TDeltaSint values is fairly small. The new, economic way to reach chemical accuracy in the calculations of the thermodynamic parameters of intermolecular interactions is proposed. In addition, interaction energy (De) and free energy change (DeltaA) for considered species have been evaluated by Carr-Parrinello molecular dynamics (CPMD) simulations and static BLYP-plane wave calculations. The free energy change along the reaction paths were determined by the thermodynamic integration/"Blue Moon Ensemble" technique. Comparison between obtained values, and available experimental and conventional ab initio results has been made. We found that the accuracy of CPMD simulations is affected by several factors, including statistical uncertainty and convergence of constrained forces (TD integration), and the nature of DFT (density functional theory) functional. The results show that CPMD technique is capable of reproducing interaction and free energy with an accuracy of 1 kcal/mol and 2-3 kcal/mol respectively.