  We show that the Franck-Condon Factor (FCF) associated to a transition
between initial and final vibrational states in two different potential energy
surfaces, having $N$ and $M$ vibrational quanta, respectively, is equivalent to
calculating the number of perfect matchings of a weighted graph with loops that
has $P = N+M$ vertices. This last quantity is the loop hafnian of the
(symmetric) adjacency matrix of the graph which can be calculated in $O(P^3
2^{P/2})$ steps. In the limit of small numbers of vibrational quanta per normal
mode our loop hafnian formula significantly improves the speed at which FCFs
can be calculated. Our results more generally apply to the calculation of the
matrix elements of a bosonic Gaussian unitary between two multimode Fock states
having $N$ and $M$ photons in total and provide a useful link between certain
calculations of quantum chemistry, quantum optics and graph theory.
