We develop a method for calculating nonadiabatic rate constants in condensed phases. The method is based on a novel path integral representation of the imaginary time flux-flux correlation function combined with an analytic continuation technique. The method is general, and may be applied to systems with arbitrarily strong coupling parameters, arbitrary anharmonic environments and any number of discrete system states. The method is applied to a nontrivial benchmark system with encouraging results. (C) 2003 American Institute of Physics.