I suppose that something crude could be attempted by looking at the basis set decomposition of the two orbitals; since they are built from the same basis set, shared bases will contribute to a Slater-type overlap according to the coefficients (or perhaps more properly the squares of the sums of corresponding coefficients) of each basis function. I am not sure how truly orthonormal practical bases in real calculations are, however, and without that strict requirement the contributions of a given basis function won't be independent of all others.
The much better method would be to run an overlap calculation in Gaussian. To my knowledge, there isn't any way to explicitly request this, but a few years ago Gaussian added the ability to calculate Franck-Condon factors and it seems to me that that is exactly what you are looking for; however, it might not give you that information without the vibrational coupling – if you try it, I would recommend paying attention to the 0-0 factor(s).
You can easily use option 11 of main function 100 of Multiwfn program to calculate overlap integral of norm of two orbitals over the whole space, its a good meaure of orbital overlap. Gaussian .fch file can be used as input file. Please check Section 3.100.11 of the manual for detail.