A multistate density functional theory in the framework of the valence bond model is described. The method is based on a block-localized density functional theory (BLDFT) for the construction of valence-bond-like diabatic electronic states and is suitable for the study of electron transfer reactions and for the representation of reactive potential energy surfaces. The method is equivalent to a valence bond theory with the treatment of the localized configurations by using density functional theory (VBDFT). In VBDFT, the electron densities and energies of the valence bond states are determined by BLDFT. A functional estimate of the off-diagonal matrix elements of the VB Hamiltonian is proposed, making use of the overlap integral between Kohn-Sham determinants and the exchange-correlation functional for the ground state substituted with the transition (exchange) density. In addition, we describe an approximate approach, in which the off-diagonal matrix element is computed by wave function theory using block-localized Kohn-Sham orbitals. The key feature is that the electron density of the adiabatic ground state is not directly computed nor used to obtain the ground-state energy; the energy is determined by diagonalization of the multistate valence bond Hamiltonian. This represents a departure from the standard single-determinant Kohn-Sham density functional theory. The multistate VBDFT method is illustrated by the bond dissociation of H2+ and a set of three nucleophilic substitution reactions in the DBH24 database. In the dissociation of H2+, the VBDFT method yields the correct asymptotic behavior as the two protons stretch to infinity, whereas approximate functionals fail badly. For the SN2 nucleophilic substitution reactions, the hybrid functional B3LYP severely underestimates the barrier heights, while the approximate two-state VBDFT method overcomes the self-interaction error, and overestimates the barrier heights. Inclusion of the ionic state in a three-state model, VBDFT(3), significantly improves the computed barrier heights, which are found to be in accord with accurate results. The BLDFT method is a versatile theory that can be used to analyze conventional DFT results to gain insight into chemical bonding properties, and it is illustrated by examining the intricate energy contributions to the ion-dipole complex stabilization.