Furthermore, not all ESDs with M S = S will contribute to the full-system states of spin S, so we can further trim the ECI basis for spin S. We accomplish this by constructing the matrix of the Ŝ2 operator in the basis of the ESDs with M S = S for each desired S separately. After diagonalizing, we obtain the excitonic configuration state functions (ECSFs)-the eigenfunctions of Ŝ2 . Then, we take only those ECSFs that correspond to the desired S, and all ESDs contributing to those ECSFs. For the details of the procedure, the notation used for ECSFs, and illustrative examples of ECSFs written in the basis of ESDs, see Appendix A. Finally, we calculate H in the spin-adapted ECI basis of ESDs, for each S separately. Before diagonalization, we transform the matrix into the basis of the ECSFs. In this way, we save some additional computational time, because the number of ECSFs is always lower than or equal to the number of ESDs for a given S, and we also directly obtain the full-system states in the most compact form, easing analysis and interpretation.