Wrong symmetry identified by cellconstructor symph package

[chenyubi14]

Dear developers,

I have encountered a weird issue. The symmetries detected by spglib and cellconstructor symph are different. (Maybe CC symph is the same as what Lorenzo meant by symmetrization with quantum espresso in discussion #76). For the identical structure, spglib identifies P6/m (175) with 12 elements, but CC symph identifies P6/mmm with 24 elements.

>>> import cellconstructor as CC, cellconstructor.Phonons
>>> import spglib
>>> 
>>> ## symmetry by spglib directly
>>> dyn = CC.Phonons.Phonons('harm_dyn', nqirr=4)
>>> symm1_spglib = spglib.get_spacegroup(dyn.structure.get_ase_atoms(), symprec = 1e-5)
>>> symm1_spglib
'P6/m (175)'
>>> 
>>> ## symmetry by CC spglib
>>> symm2_spglib = CC.symmetries.QE_Symmetry(dyn.structure)
>>> symm2_spglib.SetupFromSPGLIB()
>>> symm2_spglib.QE_nsym
12
>>> 
>>> ## symmetry by CC symph
>>> symm3_symph = CC.symmetries.QE_Symmetry(dyn.structure)
>>> symm3_symph.SetupQPoint(verbose = True)
Symmetries of the bravais lattice: 24
Symmetries of the crystal: 24
Symmetries of the small group of q: 24
>>> symm3_symph.QE_nsym
array(24, dtype=int32)

I know the correct symmetry is P6/m (175), as indicated by spglib. Currently, the code SchaMinimizer.py by default uses CC symph for symmetry identification. I am aware that spglib can be used during minimization from discussion #76 by setting use_spglib = True, but I think the current version 1.4 of SSCHA is still inevitably using CC symph.
As you can see, even though I have set minim.use_spglib = True, the execution of minim.init(verbosity=True) shows inconsistent symmetry output: 24 and 12 symmetries simultaneously.

>>> import sscha.Utilities
>>> import cellconstructor as CC, cellconstructor.Phonons
>>> 
>>> TEMPERATURE = 200 
>>> POPULATION = 1 
>>> num_conf = 50
>>> 
>>> ## load ensemble data
>>> dyn = CC.Phonons.Phonons('harm_dyn', nqirr=4)
>>> ensemble = sscha.Ensemble.Ensemble(dyn, TEMPERATURE)
>>> ensemble.load('population1', population = POPULATION, N = num_conf)
>>> minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
>>> 
>>> ## Reduce the threshold for the gradient convergence
>>> minim.meaningful_factor = 0.001
>>> minim.use_spglib = True 
>>> minim.enforce_sum_rule = True 
>>> 
>>> minim.init(verbosity=True)
Symmetry initialization...
Symmetries of the bravais lattice: 24
Symmetries of the crystal: 24
Symmetries of the small group of q: 24
Symmetry group:  P6/m (175)
Total number of symmetries: 12
Check for the imaginary frequencies...
Get the gradient for the first time...
After gradient.
>>> 

The problem lies in the SchaMinimizer.py code. Before if self.use_spglib, the code calls qe_sym.SetupQPoint(verbose = verbosity) which employs CC symph. Therefore, I cannot avoid the above symmetry problem merely by use_spglib = True.

    def init(self, verbosity = False, delete_previous_data = True, init_timer = True):
        if not self.neglect_symmetries:
            if verbosity:
                print("Symmetry initialization...")
            qe_sym = CC.symmetries.QE_Symmetry(self.dyn.structure)
            qe_sym.SetupQPoint(verbose = verbosity)

            if self.use_spglib:
                qe_sym.SetupFromSPGLIB()

                import spglib
                if verbosity:
                    print("Symmetry group: ", spglib.get_spacegroup(self.dyn.structure.get_ase_atoms()))

                self.N_symmetries = qe_sym.QE_nsym

                # Symmetrize the dynamical matrix
                self.dyn.SymmetrizeSupercell()
            else:
                fcq = np.array(self.dyn.dynmats)
                qe_sym.SymmetrizeFCQ(fcq, self.dyn.q_stars, asr = "custom")
                for iq in range(len(self.dyn.q_tot)):
                    self.dyn.dynmats[iq] = fcq[iq, :, :]

            # Save the number of symmetries
            self.N_symmetries = qe_sym.QE_nsym
        else:
            # Only apply the acoustic sum rule
            CC.symmetries.CustomASR(self.dyn.dynmats[0])
            self.N_symmetries = 1

Could you suggest a solution to this issue? Perhaps provide a new version of the function qe_sym.SetupQPoint(verbose = verbosity), to completely bypass CC symph?
Attached are my harm_dyn files for your reference.

Thanks a lot!
Yubi
harm_dyn.zip

[diegomartinez2]

I made a quick test with SrTiO3 and obtained the same symmetries with both spglib and cellconstructor symph. On the other hand with your system, cellconstructor symph doubles the symmetries obtained by spglib.
So it works for some systems and not on others?

[chenyubi14]

Dear Diego,

So it works for some systems and not on others?

I think so. In my other tests, the symmetries matched up. It is the particular attached example that shows the inconsistency.

Thank you for your prompt response. Additionally, I'd appreciate if you could take a moment to review discussion #139. I've added another question there and thought it might be overlooked.

Thanks again for your attention to these matters.
Yubi

[mesonepigreco]

I saw this problem on another system, it was due to the different threshold for symmetrization in QE and SPGLIB. In particular QE treat the system in cristalline coordinates, therefore, in case of a big cell, or big vacuum, the crystalline coordinates can be very small and much closer to higher symmetry positions. There is an open issue of some days ago. The patch is trivial, but I had not the time to publish it Il Mar 12 Dic 2023, 07:59 Diego Martinez Gutierrez ***@***.***> ha scritto:

…

I made a quick test with SrTiO3 and obtained the same symmetries with both spglib and cellconstructor symph. — Reply to this email directly, view it on GitHub <#179 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGZMNVABFA5Y7A6DMVE22VTYI76F3AVCNFSM6AAAAABAQWG3HKVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TQMRXGM3DO> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.*** com>

[chenyubi14]

Dear Lorenzo,

Thanks very much for your explanation and the forthcoming solution. Indeed, my structure is in proximity to a high-symmetry phase, which likely contributes to the issue.

I am eager to see the release of the patch. Could you please keep me updated via the GitHub discussions? Either a notification about a new version release or a direct response in this discussion thread would be greatly appreciated. Your assistance and prompt attention to this matter are highly valued.

Best,
Yubi

[mesonepigreco]

Here is the fix: SSCHAcode/CellConstructor@511f2b7

I just cherry-picked this fix into the master branch of CellConstructor.
If you are using the development version (the master branch in github), you can update CellConstructor and it should now work. Let me know (I have not tasted with your dyn).

Anyway, if you have use_spglib = True, the symmetries of quantum espresso (the symph module as you guessed is a wrapper of the quantum espresso subroutines for symmetry identification) should be ignored in the actual minimization of the code (just printed in output, could be confusing, but helps to spot potential bugs like this one. For sure, we can at least improve the clarity on which they are printed in output).

Bests,
Lorenzo

[chenyubi14]

Dear Lorenzo,

Thank you for addressing the issue with CellConstructor. The symmetry has largely improved, yet there's still a discrepancy in the bravais lattice symmetry, showing 24 instead of the expected 12. See below Symmetries of the bravais lattice: 24. Could you kindly implement another fix to this issue?

Many thanks,
Yubi

>>> import cellconstructor as CC, cellconstructor.Phonons
>>> import spglib
>>> 
>>> ## symmetry by spglib directly
>>> dyn = CC.Phonons.Phonons('harm_dyn', nqirr=4)
>>> symm1_spglib = spglib.get_spacegroup(dyn.structure.get_ase_atoms(), symprec = 1e-5)
>>> print(symm1_spglib)
P6/m (175)
>>> 
>>> ## symmetry by CC spglib
>>> symm2_spglib = CC.symmetries.QE_Symmetry(dyn.structure)
>>> symm2_spglib.SetupFromSPGLIB()
>>> print(symm2_spglib.QE_nsym)
12
>>> 
>>> ## symmetry by CC symph
>>> symm3_symph = CC.symmetries.QE_Symmetry(dyn.structure)
>>> symm3_symph.SetupQPoint(verbose = True)
Symmetries of the bravais lattice: 24
Symmetries of the crystal: 12
Symmetries of the small group of q: 12
>>> 
>>> print(symm3_symph.QE_nsym)
12

[mesonepigreco]

Hi, that is not an error, the symmetry of the bravais lattice refer to the maximum possible amount of symmetries in the primitive cell considered (ignoring the actual positions of the atoms). The real symmetries of the structure (those used by the code) is indicated as symmetries of the crystal (12 which is correct). Then also the symmetries of the small q point are written, which in the default case (q equal to gamma) corresponds to all the synmetries (12 in your case). So it is working correctly! Bests, Lorenzo Il Mar 19 Dic 2023, 23:12 chenyubi14 ***@***.***> ha scritto:

…

Dear Lorenzo, Thank you for addressing the issue with CellConstructor. The symmetry has largely improved, yet there's still a discrepancy in the bravais lattice symmetry, showing 24 instead of the expected 12. See below *Symmetries of the bravais lattice: 24*. Could you kindly implement another fix to this issue? Many thanks, Yubi >>> import cellconstructor as CC, cellconstructor.Phonons >>> import spglib >>> >>> ## symmetry by spglib directly >>> dyn = CC.Phonons.Phonons('harm_dyn', nqirr=4) >>> symm1_spglib = spglib.get_spacegroup(dyn.structure.get_ase_atoms(), symprec = 1e-5) >>> print(symm1_spglib) P6/m (175) >>> >>> ## symmetry by CC spglib >>> symm2_spglib = CC.symmetries.QE_Symmetry(dyn.structure) >>> symm2_spglib.SetupFromSPGLIB() >>> print(symm2_spglib.QE_nsym) 12 >>> >>> ## symmetry by CC symph >>> symm3_symph = CC.symmetries.QE_Symmetry(dyn.structure) >>> symm3_symph.SetupQPoint(verbose = True) Symmetries of the bravais lattice: 24 Symmetries of the crystal: 12 Symmetries of the small group of q: 12 >>> >>> print(symm3_symph.QE_nsym) 12 — Reply to this email directly, view it on GitHub <#179 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGZMNVG6MRIXATCRR32PMNLYKIGORAVCNFSM6AAAAABAQWG3HKVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM3TSMBSGM2DI> . You are receiving this because you commented.Message ID: ***@***.***>