search

SSCHAcode / python-sscha

The python implementation of the Stochastic Self-Consistent Harmonic Approximation (SSCHA).
GNU General Public License v3.0
55 stars 21 forks source link

Imaginary part error (bug in CellConstructor?) #112

Closed bkchang closed 1 year ago

bkchang commented 1 year ago

Dear Developers,

I'm trying to use your code to obtain all-real phonon dispersion starting from QE dynamical matrices with imaginary phonon frequencies. I'm following your PbTe tutorial.

I ran the following script:

dyn = CC.Phonons.Phonons("QE-dynmat/LCO.dyn", nqirr=24)
dyn.Symmetrize()
dyn.ForcePositiveDefinite()
ensemble = sscha.Ensemble.Ensemble(dyn, T0=150, supercell=dyn.GetSupercell())
ensemble.generate(N=10)

And got the following error message:

Traceback (most recent call last):
  File "/global/cscratch1/sd/bkchang/LCO-LTO/SSCHA/sscha-script.py", line 19, in <module>
    ensemble.generate(N=10)
  File "/global/homes/b/bkchang/.conda/envs/SSCHA/lib/python3.10/site-packages/sscha/Ensemble.py", line 1142, in generate
    self.init_from_structures(structures)
  File "/global/homes/b/bkchang/.conda/envs/SSCHA/lib/python3.10/site-packages/sscha/Ensemble.py", line 1069, in init_from_structures
    new_super_dyn = self.current_dyn.GenerateSupercellDyn(self.current_dyn.GetSupercell())
  File "/global/homes/b/bkchang/.conda/envs/SSCHA/lib/python3.10/site-packages/cellconstructor/Phonons.py", line 1890, in GenerateSupercellDyn
    dyn_supercell.dynmats[0] = self.GetRealSpaceFC(supercell_size, img_thr = img_thr)
  File "/global/homes/b/bkchang/.conda/envs/SSCHA/lib/python3.10/site-packages/cellconstructor/Phonons.py", line 2751, in GetRealSpaceFC
    fc = GetSupercellFCFromDyn(dynmat, np.array(self.q_tot), self.structure, super_structure, img_thr = img_thr)
  File "/global/homes/b/bkchang/.conda/envs/SSCHA/lib/python3.10/site-packages/cellconstructor/Phonons.py", line 4015, in GetSupercellFCFromDyn
    assert imag < img_thr, ASSERT_ERROR.format(imag)
AssertionError: 
    Error, the imaginary part of the real space force constant
    is not zero. IMAG=1.1554928846501804e-06

I checked the source code of cellconstructor/Phonons.py and the line causing the error is

assert imag < img_thr, ASSERT_ERROR.format(imag)

I find this quite weird as i thought the logic here should be

assert imag > img_thr, ASSERT_ERROR.format(imag)

Your kind reply is appreciated.

Best, Ben

mesonepigreco commented 1 year ago

Hi, sorry for the delay in the reply.

The value of the imaginary part should be zero, thus the assertion statement in the code is correct. The threshold is 1e-6, so 1.1554928846501804e-06 > 1e-6, and this raise the error. The reason is probably due to the fact that espresso saves dynamical matrices with single precision, and when you have a very big supercell dyn compared to the primitive cell (like the one you are using), it may give raise to this numerical problems.

I have two suggestions: 1) Consider starting with a smaller supercell (nqirr = 24 seems quite big), this also means that the calculation is way cheaper, and there are tricks to improve the convergence with the q-mesh sampling much faster than the harmonic result. 2) You can manually hack the code to increase the threshold to 1e-5 or higher, and get rid of the error.

Let me know if this is helpful, Bests, Lorenzo

bkchang commented 1 year ago

Hi Lorenzo,

Thanks for the detailed response. Yes, right after I posted this I realized you are implementing the correct logic. And indeed eventually I reduced the density of my ab initio nqirr, and am redoing the calculation now. I have other more general questions about constructing the ensemble, and will post it in discussions. Thank you so much and sorry for the confusion.

Best, Ben