Free energy landscape explorations have been performed for Cytochrome c Oxidases, aa3 from Paracoccus denitrificans and ba3 from Thermus thermophilus, interacting with small gas molecules (CO, NO, O2), as well as Xe. The calculations were carried out with thermodynamic perturbation theory, the validity of which has been examined by previous molecular dynamics calculations. This approach allows us to bypass the immense computational time required in such problems. The free energy surfaces are constructed as functions of the three Cartesian coordinates of the center of mass of the ligand and averaging over the orientation angles of the molecule. Hydrophilic/hydrophobic cavities and channels around the distal heme-a3 pocket were detected and the corresponding free energy minima and barriers were estimated. These free energy extrema permit us to extract kinetic parameters and to discuss the biochemical functions of the enzymes in conjunction with experimental results. The conserved cavities found in the two enzymes as well as in previous results of myoglobin demonstrate that topographical characteristics in the distal region of the active sites of the heme oxidase proteins are structurally stable.