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

[mesonepigreco]

Discussed in #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]

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]

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.
SSCHAcode/CellConstructor#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]

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

Error in routine frc_blk (1):
wrong total_weight

Regards,
Prasad

[mesonepigreco]

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]

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.