LO-TO splitting in SSCHA

[JonathanSchmidt1]

Hi Lorenzo,

I am working on some Insulator with SSCHA and was wondering what's the correct way to treat LO-TO splitting in SSCHA. Looking into the code during the SSCHA minimization the effective charges seem to be ignored (all the DiagonalizeSupercell() calls do not have loto set). As far as I remember from the school and tutorials as well as https://github.com/orgs/SSCHAcode/discussions/344 the charges should just be added to the dynamical matrix for the phonon calculation. However, calculating phonons afterwards with get_static_correction_along_path if I set the effective charges/dielectric constant I get "Cannot continue with SSCHA negative frequencies" and without charges it works fine but of course there is no LO-TO splitting. This is also independent of the dynamical matrix during the SSCHA calculation having charges or not and the same happens for a 3x3x3 kgrid. I am not sure if this is an interpolation issue leading to negative sscha frequencies randomly or something I am doing wrong with LOTO splitting.

What's the correct approach to deal with long range contributions in SSCHA?
In some older fortran part of the package this seems to also have been considered? (get_gaussian_weight_eff_charge)

thank you so much!
Edit
I checked also for only the k-points in .q_tot and in some cases it removed the error message in others the error message stayed.

Full output:

 
 ====================================
      Bubble static correction       
 ====================================
 
 T=   0.0 K
 k grid= 2 x 2 x 2 
 
Time for the step 1 / 8: 0.07892751693725586 s
(The preparation of the dynq: 6.103515625e-05 s)
NAT: 10
Time for the step 2 / 8: 0.07747220993041992 s
(The preparation of the dynq: 6.008148193359375e-05 s)
NAT: 10
Time for the step 3 / 8: 0.07843971252441406 s
(The preparation of the dynq: 6.127357482910156e-05 s)
NAT: 10
Time for the step 4 / 8: 0.07830023765563965 s
(The preparation of the dynq: 6.008148193359375e-05 s)
NAT: 10
Time for the step 5 / 8: 0.07835721969604492 s
(The preparation of the dynq: 6.008148193359375e-05 s)
NAT: 10
Time for the step 6 / 8: 0.07841253280639648 s
(The preparation of the dynq: 6.031990051269531e-05 s)
NAT: 10
Time for the step 7 / 8: 0.07741117477416992 s
(The preparation of the dynq: 5.888938903808594e-05 s)
NAT: 10
Time for the step 8 / 8: 0.07768917083740234 s
(The preparation of the dynq: 5.936622619628906e-05 s)
NAT: 10
Time to prepare the effective charges: 0.6253561973571777 s
Generating Real space force constant matrix...
Time to generate the real space force constant matrix: 0.01875782012939453 s
TODO: the last time could be speedup with the FFT algorithm.
/opt/conda/envs/sscha/lib/python3.10/site-packages/cellconstructor/ForceTensor.py:255: ComplexWarning: Casting complex values to real discards the imaginary part
  self.tensor[i_block, 3*na1:3*na1+3, 3*na2: 3*na2+3] = tensor[3*na1 : 3*na1+3, 3*nat2_sc : 3*nat2_sc + 3]
 
 ======================= Centering 2nd FCs ==========================
 
 Time elapsed for computing the centering: 0.014735937118530273 s
 Memory required for the centered tensor: 0.00031515955924987793 Gb
 
 ====================================================================
 
 Path printed in path_len.dat 
 
 
 dyn+odd dynamical matrices printed in sscha_plus_odd_dyn#q
 
-q-k=  [-0.00796029  0.01366599  0.00531408]     (2pi/A)
w(-q-k)=  [ -8.2079281   23.19006381  33.71361053  46.59075296  53.33268762
  74.66562548  83.52770072  90.48838458 144.94879952 145.46289886
 145.78484098 168.23288247 170.42693576 172.64895243 259.94622997
 264.39065379 267.23553533 424.08612489 425.59311172 429.42166485
 439.49489476 446.06999425 448.19673303 473.38396205 492.77026207
 494.94224868 557.2819306  560.35706874 775.14632493 839.59910076]   (cm-1)
Cannot continue with SSCHA negative frequencies
Traceback (most recent call last):
  File "/cluster/scratch/sjonathan/mirror/rpa_3/phonons_sscha.py", line 90, in <module>
    CC.Spectral.get_static_correction_along_path(dyn=dyn, tensor3=tensor3,  k_grid=k_grid, q_path=q_path,filename_st=f"v2_v2+d3static_freq_{PATH}_{args.n_points}.dat",T =0.0, print_dyn = True)
  File "/opt/conda/envs/sscha/lib/python3.10/site-packages/cellconstructor/Spectral.py", line 390, in get_static_correction_along_path
    dynq, v2_wq[iq,:] = get_static_bubble(tensor2=tensor2, tensor3=tensor3,
  File "/opt/conda/envs/sscha/lib/python3.10/site-packages/cellconstructor/Spectral.py", line 175, in get_static_bubble
    tmp_bubble = CC.Settings.GoParallel(compute_k, k_points, reduce_op = "+")
  File "/opt/conda/envs/sscha/lib/python3.10/site-packages/cellconstructor/Settings.py", line 283, in GoParallel
    result += function(x, **kwargs) 
  File "/opt/conda/envs/sscha/lib/python3.10/site-packages/cellconstructor/Spectral.py", line 142, in compute_k
    exit()
  File "/opt/conda/envs/sscha/lib/python3.10/_sitebuiltins.py", line 26, in __call__
    raise SystemExit(code)

Full code:

dyn = CC.Phonons.Phonons("final_dyn",nqirr=args.nqirr)


supercell = dyn.GetSupercell()
tensor3 = CC.ForceTensor.Tensor3(dyn.structure,
                                dyn.structure.generate_supercell(supercell),
                                supercell)
import subprocess
tensor3 = np.load(lzma.open("FC3"+postfix+".npxz",'rb'),allow_pickle=True)[()]

    
import ase
k_grid=list(args.k_grid)
PATH = args.path
N_POINTS = args.n_points

band_path = ase.dft.kpoints.bandpath(PATH,
    dyn.structure.unit_cell,
    N_POINTS)

# Get the q points in cartesian coordinates
q_path = band_path.cartesian_kpts()
x_axis, xticks, xlabels = band_path.get_linear_kpoint_axis()

CC.Spectral.get_static_correction_along_path(dyn=dyn, tensor3=tensor3,  k_grid=k_grid, q_path=q_path,filename_st=f"v2_v2+d3static_freq_{PATH}_{args.n_points}.dat",T =0.0, print_dyn = True)

[mesonepigreco]

Hi,
Sorry for my very late reply. I missed the question on the Forum.

As you mentioned, the code adds the LO-TO to the dynamical matrix before calculating the bubble correction. However, the error you get seems more related to a small imaginary frequency emerging somewhere during the interpolation.
You can try using a coarser k-grid and check if it works as expected.

It is bizarre that LO-TO interpolation gives imaginary frequencies and works fine without it.