You've probably only a found a local minimum (stability analysis only moves you off of saddle points). You can try modifying the initial guess to find the published state.

Is there a way to obtain a complex-valued initial guess density matrix?

All of the available initial guess methods in PySCF do not converge to the published state, and I suspect it may be because they all return dms of type np.float64.

I see that there is a way to perform "complex GHF" by simply changing the initial guess to dtype=np.complex128 (and maybe adding some small imaginary perturbations), however this still doesn't seem enough to converge onto the published state.

Sorry, I don't have enough experience but maybe someone else can comment.

You may achieve this uniformly or randomly add a small imaginary part to a current real valued densitry matrix.

With complex initial guess (like @xubwa suggested) I can get a few complex solutions, with non-zero mx and my population on atoms, although still not equals the value in reference.

However, the ghf_stability function may not be able to detect complex GHF internal instability or real2complex instability, I believe. So it still remains a question how we can tell the complex solution is right.

@jeanwsr Are your solutions noncollinear? Meaning, is there no way to make the mx and my components vanish by choosing a different set of axes? For instance, if you rotate the OP's two local-minima plots 120 degrees, you'd have two equivalent local minima that would appear, at first glance, to be properly "complicated" (mx and my and imaginary components of the density matrix). But obviously those two are still stuck in the same high-symmetry space that we are trying to escape here to reach the true global minimum.

This type of spin frustration may be something that just ought to be an initial guess option for GHF (i.e., make all spin vectors radiate away from or point towards a given coordinate).

Thanks @MatthewRHermes . I checked my previous tests and found most of the complex solutions are collinear. After managing initial guess very carefully I can get the "radiating away" spin vectors, although I'm not sure my population function is right.

Notebook attached
cghf.ipynb.zip

I think that's right. The population function used in OP's reference paper is more complicated, but I don't think it matters. It looks like your magnetization is C2v, rather than D3h, which I thought was a problem at first, until I realized that your geometry is C2v. Check out my version, which contains a density-matrix symmetry-breaker that I think would be generally useful (although the guess density matrices are still coplanar).

So I'm just learning this myself, but it turns out that since the true optimum GHF density is part of a huge degenerate manifold (see DOI:10.1063/1.4918561) in which spin vectors don't seem to point away from or towards anything in particular, but in fact have exactly the same energy as the "correct" density matrix I posted above. With that in mind, it turns out that the default SCF guess followed by one round of ghf_stability as currently implemented in PySCF is actually sufficient to get to the global minimum. You just can't see it because of the massive degeneracy.

Thanks everybody for your help with this.
I've attached my workbook below (it includes the Hirshfeld population scheme similar to the original paper) with my initial calcs + some of the suggestions above. Although nothing perfectly matches the published figure, it seems that the states with spin vectors angled 120 deg from eachother may be degenerate, in line with @MatthewRHermes last comment.

li_trimer_ghf.zip

Looks like ghf_stability is still needed to give the lowest energy in @mc-rohan 's results. However, the current ghf_instability implementation may not be correct in case of complex ghf. It should diagonalize ((A B) (B* A*)) as DOI:10.1063/1.4918561, 10.1063/1.434318 stated but currently diagonalize A+B.

Dang, I think you're right. I just noticed that there's a separate real2complex Hessian function for UHF external stability in pyscf.scf.stability. But this also implies that Newton GHF does not actually work!!