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

error while interpolation of 2x2x2 SSCHA dyn_mat using a 4x4x4 harmonic dyn_mat #79

Open mesonepigreco opened 2 years ago

mesonepigreco commented 2 years ago

Discussed in https://github.com/SSCHAcode/python-sscha/discussions/74

Originally posted by **matukumilli** March 29, 2022 Hi, The workaround suggested by Dr. Lorenzo in an earlier thread worked well. In order to avoid SSCHA -ve frequencies at finite q points one has to interpolate SSCHA dyn_mat using a harmonic dyn_mat calculated on a finer q-mesh than SSCHA q-mesh. However, when I try this interpolation way to avoid -ve freq error, for tetragonal SrTiO3, I getting this below error: .... ... Iter # 128 ==== Sum on 3rd= 0.844989E-12 Imp. ASR on 3rd: => delta FC= 0.844988E-12 Imp. permut sym: => delta FC= 0.800810E-12 Convergence reached within threshold: 0.100000E-11 Total FC relative variation: 0.128277E+00 Time elapsed for imposing ASR: 186.2998547554016 s ============================================================ WORKING ON: [0. 0. 0.] SUPERCELL SIZE: 2 2 2 NAT: 10 The unit cell: 5.5080868904697313 7.3440103850181005E-003 0.0000000000000000 3.0484109872987837E-002 5.5080073591324155 0.0000000000000000 -2.7692854120418908 -2.7576755993832940 3.8814121235103820 Total weight: 7.6666666666666670 NR1,2,3 = 2 2 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Error in routine frc_blk (1): wrong total_weight %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% stopping ... STOP 1 I couldn't able to figure out the reason for this. Need your help. Original harmonic dyn_mat calculated on 2x2x2 q-mesh; to avoid the -ve freq in SSCHA dyn_mat, I tried to interpolate the difference between SSCHA and harmonic dyn_mat using harmonic dyn_mat on 4x4x4 q-mesh. *************** Input script: *************** from __future__ import print_function from __future__ import division import cellconstructor as CC import cellconstructor.Phonons import sscha, sscha.Ensemble import cellconstructor as CC import cellconstructor.ForceTensor import cellconstructor.Structure import cellconstructor.Spectral import numpy as np import ase from ase.calculators.espresso import Espresso from ase.visualize import view dyn = CC.Phonons.Phonons("../dyn_pop3_",6) supercell = dyn.GetSupercell() tensor3 = CC.ForceTensor.Tensor3(dyn.structure, dyn.structure.generate_supercell(supercell), supercell) d3 = np.load("d3_realspace_sym.npy")*2.0 tensor3.SetupFromTensor(d3) tensor3.Center(Far=2) tensor3.Apply_ASR() k_grid=[10, 10, 10] path = ase.dft.kpoints.bandpath("GXPNGZS1", dyn.structure.unit_cell, npoints = 150) qpath = path.cartesian_kpts() harm_dyn = CC.Phonons.Phonons("../../../matdyn", 6) harm_dyn_fine_grid = CC.Phonons.Phonons("matdyn_4x/matdyn", 24) new_cell = harm_dyn_fine_grid.GetSupercell() big_dyn = dyn.Interpolate( dyn.GetSupercell(), new_cell, support_dyn_coarse = harm_dyn, support_dyn_fine = harm_dyn_fine_grid) big_dyn.ForcePositiveDefinite() CC.Spectral.get_static_correction_along_path(dyn=big_dyn, tensor3=tensor3, k_grid=k_grid, q_path = qpath, filename_st="v2_v2+d3static_freq-GXPNGZS1.dat", T =50.0, print_dyn = False) Regards, Prasad
matukumilli commented 2 years ago

Hi Lorenzo, due to large size of d3 tensor file, I failed to attach here, rather I sent them in email. Pls use those files to reproduce the error message.

Regards, Prasad

mesonepigreco commented 2 years ago

The issue seems not related to the d3 part, but to the interpolation using support dynamical matrices. I created a new branch of cellconstructor which fixes this error. https://github.com/SSCHAcode/CellConstructor/pull/48

However, the code seems not to be able to reproduce the same output as the previous interpolation and the testsuite fails (up to today), therefore more testing is needed.

matukumilli commented 2 years ago

Similar behaviour is observed at little lower temperature (T=20K) as well:

Error in routine frc_blk (1): wrong total_weight

Regards, Prasad

mesonepigreco commented 2 years ago

Thanks for the report. I'm working on it. It seems a problem with the interpolation using a support matrix. We will fix this soon, however, if you want to proceed with the calculation, you can try to interpolate without the support matrix and use the ForcePositiveDefinite after the interpolation and before passing it to the get_static_correction_along_path (as you are doing), it should avoid errors in the bubble calculation due to imaginary frequencies.

matukumilli commented 2 years ago

dyn.ForcePositiveDefinite() before passing it to the get_static_correction_along_path() didn't help: q= [ 2.76026848e-03 2.76026848e-03 -1.06370064e-19] (2pi/A) w(q)= [ -4.03429147 -3.46423257 8.83285856 44.34200573 47.2531449 110.37484302 111.83469055 136.15802509 140.4898301 140.75631651 146.07548221 171.62276582 173.40597864 175.80646154 249.52104453 249.74820731 261.18240083 425.39532056 427.77205168 428.40073337 429.5592279 432.71955735 433.20384253 455.8704559 488.9355318 497.85667105 532.92236704 532.94816711 772.08949798 825.88642211] (cm-1) Cannot continue with SSCHA negative frequencies

I thought by doing so it only fixes at high symmetry points.