Ab initio calculations are reported of the rate of tautomerization by double-hydrogen transfer of porphine and three of its isotopomers. Both synchronous tone-step) and asynchronous (two-step) hydrogen tunneling mechanisms are considered. Geometries and force fields are calculated at the stationary points by means of a nonlocal density functional method that yields accurate equilibrium structures and vibrational spectra. Potential-energy surfaces are constructed in terms of all 73 in-plane normal-mode coordinates at the transition state, the mode with imaginary frequency being taken as the reaction coordinate. Hydrogen tunneling calculations are performed by means of a simplified instanton method that has proved reliable in calculations on smaller systems. The full multidimensional potential is used, and adiabatic separation of the normal modes from the reaction coordinate is avoided. The coordinates of the transverse modes are coupled linearly to the reaction coordinate and all modes are allowed to mix freely with each other along the reaction path. Direct evaluation of the instanton path is not necessary. To calculate the tunneling rate constant, it is sufficient to evaluate the one-dimensional instanton action along the reaction coordinate and to correct it for coupling with transverse vibrations. This makes the method computationally very efficient compared to other multidimensional approaches. For the synchronous mechanism, the calculations closely follow the previously established procedure, but for the asynchronous mechanism, generalization to an asymmetric barrier is required. This is achieved by dividing the normal-mode displacements that determine the couplings into symmetric and antisymmetric components which enhance and suppress the tunneling rate, respectively. The relative energies at the stationary points of the density-functional potential are calculated both by density functional theory (DFT) and by the Hartree-Fock method at the DFT geometry. The two methods yield results that are quite different. Comparison with a large set of experimental data comprising four isotopomers and a wide range of temperatures, indicates that neither method yields accurate energies but that some adjustment of the barrier height and the cis-trans energy difference is necessary to obtain satisfactory rate constants for the asynchronous mechanism. The other calculated parameters are used without adjustment. All parameters are combined to construct the potential required for the instanton calculations. A good fit to all available kinetic data is obtained, indicating that the method accounts accurately both for the isotope and the temperature dependence of the rate of tautomerization, It is shown that, in order to achieve this result, it is essential to include all linear couplings, since the balance between symmetric couplings, which enhance the tunneling rate, and antisymmetric couplings, which suppress it, varies between isotopomers. All dynamics calculations are performed with a newly developed code, which is designed to use the output of standard quantum-chemical codes and requires only minutes of CPU time on a standard workstation.