While PCM is useful, it cannot explicitly describe the solvent structure and thermal fluctuations, leading to inaccuracies in modeling intermolecular interactions. This is clearly demonstrated by the calculation with explicit solvent molecules (Table ). To address this, the RISM method or the recently developed fluctuating charge force field, available in OpenMolcas, could prove beneficial. Furthermore, the current formulation does not yield the variational minimum of the state-averaged energy. Given these limitations, it would be worth exploring the possibility of integrating the recent formulation for RISM. Additionally, the present method cannot be applied to molecular dynamics simulations that involve non-adiabatic transitions, as it assumes the solvent is polarized with respect to a specific density. A preliminary implementation of a dynamic solvation model, including analytic derivatives, seems promising and could be a potential direction for future work. We plan to formally release the developed code as a part of OpenMolcas in the near future.