Suppose we have a parametrized molecular Hamiltonian $H(R)$, where $R$ are the nuclear coordinates. If we let $B(R)$ be some chosen atomic basis set (like STO-3G) centred at $R$, then using $B(R)$, we let $C(R)$ denote the transition matrix from atomic to molecular orbitals at $R$.
In general, there are many different choices for $C(R)$ at a given $R$, as HF can converge to different unique solutions. However, I want $C(R)$ to vary continuously over some range of $R$. In general, is it possible to enforce this condition when performing Hartree-Fock? If so, how?
I tried starting at a set of parameters $R_0$, computed $C(R_0)$, and then computed $C(R_0 + \Delta)$, $C(R_0 + 2\Delta), ...$ iteratively by using the density matrix yielded by the previous Hartree-Fock step, but this didn't seem to work: I'm still getting some large jump discontinuities in $C(R)$ at certain values of $R$.