Controlling ⟨Ŝ2⟩ in broken-symmetry density functional theory calculations via constrained optimization

Accurate determination of magnetic exchange coupling constants (J) from density functional theory (DFT) remains challenging, particularly for open-shell systems where broken-symmetry (BS) solutions suffer from spurious spin contamination that systematically exaggerates J values. Several methods have been proposed to address this problem by adjusting the mapping scheme from the DFT energies to the Heisenberg-Dirac-van Vleck effective spin Hamiltonian energies. In this work, we explore a different route by imposing a constraint on the DFT energy that enforces a target value of the spin-squared expectation, ⟨Ŝ2⟩, using a Lagrange multiplier approach. By explicitly controlling the spin character of the electronic state, the method attempts to overcome limitations of standard BS calculations to describe magnetic interactions. As part of the theoretical formulation, we derive analytical expressions for the gradient of the spin-squared expectation value with respect to the spin-resolved density matrices, which are required for the practical implementation of the constraint within a generalized Kohn-Sham scheme. These expressions are general to any single-determinant method and remain valid for arbitrary spin states. We apply the spin-constrained approach to the calculation of J couplings and compare with three energy-difference-based schemes for a set of representative systems, including H2He, H3He3 arranged in an equilateral triangle, and a bis(μ-hydroxo) Cu(II) complex. Across all cases, the constrained formulation yields systematically lower and more consistent exchange couplings across different density functional approximations. This work establishes a robust and general route for incorporating spin-state constraints into DFT-based studies of magnetic exchange interactions.