We develop an algorithmic, metaheuristic approach to the definition of molecule-fixed axes orientation in molecules of arbitrary size. The goal is to simplify the treatment of overall rotation and rotation–vibration interaction in rovibrational Hamiltonians. Considering the kinetic elements of the rovibrational Hamiltonian, given by the G matrix, we select the optimal orientation of molecule-fixed axes minimizing specific G matrix elements. To such an end, we develop a global minimization method based in a hybrid Simulated Annealing+Powell's local minimization. The parameters of the method are calibrated using a set of non-rigid molecules: Acetaldehyde, glycolaldehyde, methyl formate and ethyl methyl ether. The results show that the principal axes of inertia do not give the simplest form to the pure rotational contribution. However, minimization of the G matrix rotational element does. Finally, we observe that in the cases considered it is not possible to nullify all the rotation–vibration coupling elements, since the torsional motions are coupled with the overall rotation. However, the treatment yields the optimal solution. The methodology proposed allows also to simplify simultaneously the pure rotational+rotation–vibration coupling elements in rovibrational Hamiltonians.