Could you post some example which shows breaksym=True is problematic?

Could you post some example which shows breaksym=True is problematic?

Here is an example where RHF converges normally, but GHF fails to converge.

basis = "def2-svp"
atom = """
C                 -0.77452336   -0.08132633    1.74172506
C                 -0.26411504   -0.11129929    0.43692360
C                 -0.35397350    1.02601929   -0.37692764
C                 -0.95424489    2.19330919    0.11402082
C                 -1.46465579    2.22328125    1.41882131
C                 -1.37479443    1.08596370    2.23267366
H                 -0.70591476   -0.94969444    2.36311840
H                  0.03573621    1.00313485   -1.37317189
H                 -1.02285582    3.06167644   -0.50737346
H                 -1.92297618    3.11453305    1.79367143
H                 -1.76450608    1.10884747    3.22891718
O                  0.34840202   -1.30241314   -0.06404605
C                  1.35270587   -0.95567870   -1.02115052
H                  2.09459791   -0.34266532   -0.55348523
H                  0.90396901   -0.41799437   -1.83011950
C                  2.01233961   -2.23841673   -1.56065630
O                  2.09302887   -3.25968240   -0.82985109
N                  2.54773286   -2.26985691   -2.92932924
H                  3.50527213   -1.98172537   -2.91939174
C                  2.45289144   -3.63693794   -3.46130165
C                  3.06213638   -3.98156820   -4.81429695
C                  3.20775770   -4.83278916   -2.92168779
H                  1.41330079   -3.88850264   -3.43186598
H                  2.21258260   -3.95040754   -5.46405646
O                  3.42944279   -5.30612203   -1.77696850
N                  3.67522246   -5.14821806   -4.31096996
C                  3.66579173   -6.10484804   -5.36714016
H                  4.49655524   -6.75828572   -5.20055325
C                  3.84517881   -5.06783193   -6.57540908
C                  4.72191908   -5.64179240   -7.70390186
H                  4.27386857   -6.53609107   -8.08387642
H                  4.80563701   -4.92235784   -8.49149645
H                  5.69541534   -5.86571901   -7.32041536
C                  2.45462245   -4.72770184   -7.14309246
H                  2.00895589   -5.61186346   -7.54872716
H                  1.83452668   -4.34238508   -6.36084427
H                  2.55421902   -3.99253283   -7.91413529
S                  4.42031248   -3.49330876   -5.80985222
C                  2.41628932   -6.97906376   -5.58182741
O                  1.24772273   -6.54739905   -5.20566321
O                  2.52902641   -8.14944067   -6.13940253
H                  3.39110571   -8.46788950   -6.41690875
"""

from pyscf import gto
from pyscf import dft
from pyscf import scf



mol = gto.M(atom = atom)
mol.max_memory = 100000 # Allocate 100GB memory
mol.basis = basis
mol.build()
print(f"mol.nbas = {mol.nbas}")
print(f"mol.nao = {mol.nao}")
print(f"mol.natm = {mol.natm}")
print("\n")

mf = scf.GHF(mol)
mf.max_memory = 100000 # Allocate 100GB memory
mf.verbose = 4
energy = mf.kernel()

mf = scf.RHF(mol)
mf.max_memory = 100000 # Allocate 100GB memory
mf.verbose = 4
energy = mf.kernel()

Your change is breaking quite a few tests. Are the results from new guess more correct or not correct in these cases?

The brokensym initial guess for GHF is quite tricky when brokensym is really needed. There're some discussions in #3058, #2591, etc. And your issue is another problem that brokensym may be a trouble when brokensym is not needed.

In the "no-brokensym" condition I think using mol.RHF().run().to_ghf().run(), or something similar with UHF, can avoid the GHF initial guess problem. And this makes sense since it indicates that the user has already known there's no brokensym, while mol.GHF() sounds like the user does not know whether there's brokensym and let it decided by program.

However, it still remains a question whether the default behavior of mol.GHF() should be brokensym or not. Setting no brokensym as default may cause troubles when intrinsic brokensym exists.
I think it depends on whether the new guess can pass the tests (or fix the tests if previous tests are wrong).

Your change is breaking quite a few tests. Are the results from new guess more correct or not correct in these cases?

The brokensym initial guess for GHF is quite tricky when brokensym is really needed. There're some discussions in #3058, #2591, etc. And your issue is another problem that brokensym may be a trouble when brokensym is not needed.

In the "no-brokensym" condition I think using mol.RHF().run().to_ghf().run(), or something similar with UHF, can avoid the GHF initial guess problem. And this makes sense since it indicates that the user has already known there's no brokensym, while mol.GHF() sounds like the user does not know whether there's brokensym and let it decided by program.

However, it still remains a question whether the default behavior of mol.GHF() should be brokensym or not. Setting no brokensym as default may cause troubles when intrinsic brokensym exists. I think it depends on whether the new guess can pass the tests (or fix the tests if previous tests are wrong).

Thanks for the feedback.

You are right, I noticed that this commit broke a number of unit tests. I am currently investigating the root causes of these failures. Apart from the tests that verify the initial guess itself, some errors seem to stem from deviations in the subsequent calculations caused by the change in the guess. This certainly highlights the sensitivity of the results to the initial guess.

I also agree with your point that hard-coding brokensym=False as the default might be too arbitrary, as it could suppress intrinsic symmetry breaking when it is actually needed. Instead of a hard switch, I think exposing an interface or argument (e.g., passing a parameter to control the symmetry breaking behavior) might be a more reasonable approach.

I will update the PR once I verify whether the currently "broken" tests are producing physically incorrect results or just different ones.