In this paper we develop a “diatomic in molecules semiempirical ligand field” (DIMSELF) method to calculate ground and excited many-body potential energy surfaces for an arbitrary transition metal ion in an arbitrary complex system. This method is not restricted to a high-symmetry environment and is meant to be inexpensive and suitable for nonadiabatic excited states dynamics on-the-fly. Within the approximations employed, the method includes full CI (configuration interaction) and SO (spin-orbit) interactions, essential to the description of nonradiative transitions such as those of myoglobin in the presence of carbon monoxide. We test our method against high level ab initio calculations for a simple model system of myoglobin’s heme pocket. Finally, we discuss our results and compare with previous calculations in the literature. Transition metal complexes play a fundamental role in many technological and biological areas. They act as catalysts in the oil industry during the processes of cracking and reforming, they are also fundamental for oxygen transport in blood, and they are key to plant respiration. Some of these processes can be well understood without the need to include the small contribution of spin-orbit (SO) interaction which, for the most part, is of the order of the error one performs in ab initio calculations of these kind of systems. In some cases, however, nonradiative transitions occur between states with different spins. Of particular interest are heme proteins, a biologically important group of molecules that have a unique common active site: an iron-protoporphyrin-IX complex, where changes in spin state can be an integral part of the protein function. One of the most striking biological manifestations of this modulation is the binding of molecular oxygen and carbon monoxide (CO) 1-9 to the five-coordinated ferrous heme site of globins. Another wellstudied example is the situation of several spin changes in the enzymatic activation cycles of cytochromes P450. These transitions are solely due to SO coupling and are otherwise forbidden. Ab initio calculations are performed within the framework of Hund’s Case A and hence assume that S, the total spin of the system, is a good quantum number. Spin-orbit coupling matrix elements could then be added, and rediagonalization of the Hamiltonian matrix would provide us with the spin coupled eigenstates of the system. Obtaining several electronic excited states for a given spin multiplicity is complex, and for mediumsized systems only configuration interaction for single excitations (CIS) is affordable. However, a more complete description of excited states might be required for the study of dynamics. For large systems, where the environment can be introduced by mixing quantum chemistry with molecular mechanics, electronically excited states are rarely computationally accessible and only the fundamental ground state for each spin multiplicity might be computed. For nonsymmetric systems, even at the density functional theory (DFT) level of theory, several attempts with different initial guesses should be considered in order to achieve the true ground state. The computational effort is drastically increased when several nuclear configurations are required for a complete study of the dynamics, and calculations on-the-fly, including excited states and SO coupling, are simply out of the question. Therefore, a semiempirical approach that includes excited states and SO coupling is very desirable. The method we develop in this paper is a generalization of Gerloch’s 10 cellular AOM (angular overlap model) version of ligand field theory, combined with a semiempirical version of Ellison’s 11-13 DIM (diatomics in molecules). Neither ligand field theory nor the semiempirical version of DIM are exact ab initio methods; not all the parameters can be obtained by analytical methods and they ought to be fitted either to experimental data or other calculations. On the other hand, once these parameters are known, fast computation of ground- and excited-state potentials can be performed without having to recalculate wave functions and their matrix elements. We test this method against a simplified model of the heme