The external potential with the contribution of the MM environment is plugged into the KS-DFT Hamiltonian, and the eigenvalue problem is solved to determine the optimal QM density. Here, the complexity arises due to the induced electrostatic potential v ind MM since it is computed through the induced MM dipoles, which depend on the QM density. This means that the induced potential is not stationary with respect to QM density variations. Thus, for a fully-polarizable KS-DFT/AMOEBA model, solving for the induced dipoles at each variation of QM density is required until self-consistency.