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

[matukumilli]

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

[mesonepigreco]

Ho, sorry for the delay in my reply.

It seems there is an error when applying the acoustic sum rule. The code does 128 iterations, which is really too much.
This is probably due to a wrong centering of the. Matrix.
Try to increase the Far = 3 (or higher) when centering the third order force constant.

tensor3.Center(Far=3)

[matukumilli]

Hi Dr. Lorenzo,

Thanks for your reply. I tried with Far =3 and 5, the same error still persists. Also no effect on the number of iterations to converge the ASR on 3rdIFCs. This slow convergence is also noticed in the cubic phase of SrTiO3 but no issue of SSCHA -ve frequencies, so I didn't worried about the interpolation using finer harmonic phonon q-grid. However, in tetragonal SrTiO3 at 50K, I am facing this SSCHA -ve freq issue despite a good free energy minimisation.

Any further help on this would be highly appreciated.

Regards,
Prasad

[mesonepigreco]

Hi, what happens if you comment the tensor3.Center(Far = 2) line?
(the interpolation will not produce meaningful results, but lets spot if the error is on the centering)

[matukumilli]

The error with commenting the line : #tensor3.Center(Far=5) is as follows:

Iterative ASR imposition with pow= 0.000

Initialize indices for permutation...
ERROR: new triplets found during the permutation symmetry initialization
the execution stops here.
If the 3rdFC is centered, try to repeat the centering with higher Far
If the 3rdFC is not centered, try to repeat the ASR imposition with PBC=True

[mesonepigreco]

I open an issue related to this. Sorry for the inconvenience, we will try to workout what is happening.

[zhaowb12138]

I got this error too

think the error is from PHonon/PH/SUBROUTINE frc_blk
If you use the 222 dynamic matrices perform the operation of drawing the phonon spectrum of QE
[q2r.x <q2r.in>q2r.out
matdyn.x < matdyn.in > matdyn.out]
you will get the same error
The solution is change [/home/zwb/qe-6.8/PW/src/wsweight.f90/function wsweight(r,rws,nrws) parameter (eps=1.0d-6)] from 1.0d-6 to 1.0d-7
This is the code where the error occurs

In my test,I found in one point the right way is
ck=-5.403862182618013E-006 < 1e-6 not return
abs(ck)=5.403862182618013E-006 > 1e-6 nreq=nreq
but in my dynamic matrix
ck = -4.229048151671577E-007 < 1e-6 not return
abs(ck)= 4.229048151671577E-007 < 1e-6 nreq=nreq+1
So the weight is wrong.
I think changing the parameter eps from 1e-6 to 1e-7 will not the correctness of the program. what do you think? And how can I change the parameter eps in SSCHA? It seems to be packaged in a symph.
And there is another question when I test, I use two sets of dynamic matrices. The first group is the dynamic matrix for the convergence of SSCHA dyn_pop5_. Another group is calculated from QE matdyn , and uses the same structure from dyn_pop5_(I just copy the structure message from dyn_pop5_ and put them into pw.x ).
Then something happened, the atomic positions of the dynamic matrix calculated by QE is different from dyn_pop5_*

The difference is not very big, But this difference lead to matdyn* calculate normally,but dyn_pop5_* not.
The same thing happened in dyn_pop5_* and the OUT.put file

I feel like such a small difference can be counted as a numerical error, but since it is the cause of this error, I think maybe I should ask you ,just for check

[mesonepigreco]

Hi, thanks a lot for this report.

the wsweight function is inside the FModules/from_matdyn.f90 file (the directory FModules contains the source compiled into the symph module), I just committed the change of the epsilon you proposed in the master. You are right, this value should be as small as possible. I will investigate this error in the structure, it could be a problem on the conversion between angstrom to Bohr, for which the Ry constant implemented in quantum espresso could be a bit different than the one used in the sscha code. I remember that ASE, when reading espresso output, initializes the Ry constant to a conversion slightly different than the standard one.