Interpolation of Hessian

[tveskaeg]

Dear developers,
I've gone through your tutorials and I'm currently trying to do a bit more on H3S (200GPa, Im-3m phase with 4 atoms/uc). So far, I've done 5 populations with 100, 100, 200, 500, and 1000 individuals in the 2x2x2 supercell (blue, red, yellow, purple, green in the plots). The other parameters follow loosely the tutorial.

These are my convergence plots:

After each population I calculated the Hessian (3rd and 4th order). The following shows phonon dispersions interpolated with QE in the Hessian of the two last populations (3rd and 4th order are almost the same, harmonic starting point as dashed line):

Now, I'm trying to interpolate the Hessians as described in the PbTe tutorial, but I encounter some problems there: In the first check with same k/q grid, which should return the same Hessians, I find these differences in the frequencies:

And further, the interpolated Hessian/dynq files don't contain the part for the symmetrically equivalent q points.

So, my questions:

• Why does the interpolation not return the original Hessians? Are the differences in frequency to be expected or are they maybe due to a not fully converged calculation?
• Is it possible to get QE dynq files that contain all the equivalent points (without running the interpolation routine for each point manually and pasting together the files)? Unfortunately, the QE post-processing needs all points explicitly.

Thanks and all the best!

P.S.: I updated my SSCHA and CC versions just this week, so I should have the latest release.
P.P.S.: This is my interpolation script: hessian_interpolate.py.txt

[mesonepigreco]

Hi,
Thanks for reporting this detailed test.

I tag @rafbianco in the discussion, as he is the main developer of the interpolation part, and I believe he can have more insight than me on why the interpolation is not giving exactly the same result as the hessian method (you need to compare with the one without the 4th order, as up to now we can only interpolate the 3rd order)

I can immediately answer why the code is not generating the full dynamical matrices. The reason is that you are computing the interpolation along a path, so the code is explicitly computing and saving only the q points you have expressed (if you use a longer path like the one in which you plotted the dispersion, it would occupy a lot of memory to save all the symmetry equivalent q points).

There are two workarounds for that:

1. You can use the quantum espresso q2qstar.x utility (from PHonon), which reads only the first dynamical matrix and recomputes the full dynamical matrix of the star of q.
2. You can use the subroutine:

CC.Spectral.get_static_correction_interpolated(dyn, tensor3, T, new_supercell, k_grid)

Instead of get_static_correction_along_path, where the new_supercell and k_grid are (2,2,2). new_supercell defines the q mesh of the output dynamical matrix, k_grid is the interpolation used to integrate the bubble diagram of the self-energy.

This will automatically return a CC.Phonons.Phonons() object that you can save with save_qe("hessian"), and it will generate the full dynamical matrix, with all the q points commensurated.
However, I think there is a small bug in the subroutine currently in the master (already fixed in the fast_upsilon branch that will be merged very soon) of CellConstructor, so this second option works at the time I'm writing only with the fast_upsilon branch.

I hope to be of help,
Bests,
Lorenzo

[tveskaeg]

Dear Lorenzo (@mesonepigreco),

as you suggested interpolating only the SSCHA/harmonic differences for a better starting point, I'm coming back to you with a similar topic: In the 2016 Nature article on which the H3S tutorial is based, there's a computational detail stating that the difference between harmonic/anharmonic dynamical matrices was interpolated from 3x3x3 to 6x6x6 and added to the harmonic 6x6x6 matrices to perform the electron-phonon calculation (-> isotropic Migdal-Eliashberg, so I assume QE/EPW?).
For further calculations one should use the physical matrices (the Hessians), so which difference could I interpolate in that case? Is it maybe performing the Hessian interpolation on files obtained with get_full_hessian=False and adding them to the corresponding harmonic matrices?

Best,
Roman

[mesonepigreco]

Dear Roman,
In that case, the best thing to do is to interpolate the SSCHA auxiliary dynamical matrix using the harmonic ones as you explained on a bigger cell, then use the auxiliary dynamical matrix you get to compute the physical phonons (the hessian matrix) on the desired supercell with the Spectral module, employing the third-order force constant file (d3_sym.npy) in the supercell you computed with get_free_energy_hessian() and a big k-mesh.
This is the most precise way of doing things.

Indeed, an easier approach is interpolating the Hessian dynamical matrices directly, adding and subtracting either the harmonic or the SSCHA auxiliary dynamical matrices in a bigger mesh you get from get_free_energy_hessian(), but it is less precise, and requires bigger starting supercells to converge.

Bests,
Lorenzo

[tveskaeg]

Dear Lorenzo,
I'm testing the procedures you have suggested, but I'm encountering some problems.
In the included python script, I want to interpolate the SSCHA dynamical matrices from a 2x2x2 to a 3x3x3 grid (it's all based on the 5th population in the H3S calculation I've shown here). I'm trying both ways, the dyn.Interpolate routine with and without harmonic support dynamical matrices.
SSCHA_222_to_333_interpolate_py.txt
(dyn_end_population5* are the final SSCHA matrices, matdyn_333_* and matdyn_222_* are the harmonic QE files for the corresponding grid)

1. without: this routine seems to work quite well, but the order of the resulting q-points is not the QE order, and for higher grids (like 6x6x6) it produces some other points than QE does in that grid, which do not even appear in the star (although they could be somehow symmetrically related or maybe it's just a printing issue -> manually reordering and using the new points' files with the QE dynq0 gives me the same dispersion).

2. with support: the resulting files are basically just the support harmonic matrices on the fine grid. It seems like that the calculated FC difference to be added is just (almost) zero. However, the q order is correct here (since it's taken from the support files, as far as I saw).

Could you maybe just quickly check my script (it should be a min. WE), if you can spot an error or if there's something with the routines?

Concerning the Hessian, so I should calculate the d3_sym.npy and a tensor3 object from the last ensemble in the coarse supercell, then interpolate the last SSCHA matrices (either with or without harmonic support) on a fine grid, and finally construct a dyn = CC.Phonons.Phonons object with those, on which I use CC.Spectral.get_static_bubble / CC.Spectral.get_static_correction, right? Or is there another routine?

All the best,
Roman

[tveskaeg]

Dear Lorenzo (@mesonepigreco),
I was just wondering if you saw my last post (the one above this: #49 (reply in thread))

I don't think I understood correctly how to interpolate the Hessians "in the most precise way". Do I need to construct a tensor3 object in order to use the Spectral module?
I'm struggling explicitly with the size of 3rd-order force constants. I create the variable d3 = np.load("d3_realspace_sym.npy")*2.0 from the file "d3_realspace_sym.npy", which is written during the call of ens.get_free_energy_hessian(...), where ens is my final ensemble in the coarse supercell. Hence, the dimensions of d3 are (n_mode, n_mode, n_mode), where n_mode is the number of modes in the coarse supercell (for H3S, 4 atoms in the uc, 2x2x2 grid: 3 * 4 * 2^3 = 96).
If I setup the tensor3 object for the 2x2x2 supercell, I can feed the d3 via tensor.SetupFromTensor(d3).
If I try it on the higher supercell, I'm obviously missing some data (as n_modes for a 3x3x3 cell would be 3 * 4 * 3^3 = 324).

And in the end, which method from the Spectral module should I use here?

All the best,
Roman

[mesonepigreco]

Hi Roman,
Sorry, As I marked as answered, I did not spot your comment. Sorry for the delay in my answer.
When you use Interpolate with a support matrix, you must use the harmonic dynamical matrix as support, not the SSCHA minimized matrix.
In your script, the following line has a problem:

new_dyn_support = dyn.Interpolate(coarse_grid=dyn.GetSupercell(), fine_grid=support_dyn_fine.GetSupercell(), support_dyn_coarse=dyn, support_dyn_fine=support_dyn_fine)

In particular the

support_dyn_coarse=dyn

That is the problem. In this case, the code is removing from dyn to dyn, and you interpolate zero. You should replace support_dyn_coarse with the harmonic dynamical matrix in the same q mesh as the dyn Phonon class, which in your case is support_dyn_coarse. You should replace the line as:

support_dyn_coarse=support_dyn_coarse

And it should work.

Regarding the third-order tensor, indeed, the supercell in which you define your Tensor3 object must match the d3 calculations. However, since the d3 is then stored in real space after loaded, when you use the Center() method, it will no more belong on the original supercell (it is in an infinite cell where all the elements whose perimeter between the three atoms that identify each element of the d3 not included in the original supercell are 0).

Thanks to this real space parameterization of Tensor3, we can get the tensor at any q point by computing the Fourier transform (even q points non-commensurate with the original supercell).

[rafbianco]

Dear tveskaeg,

thank you for the analysis. Lorenzo has already replied to many of your questions. Regarding the specific point of the mismatch, this should be due to the ASR imposition. This operation inevitably changes the FCs. To double check that, do the test without imposing the ASR. As long as you use the grid commensurate with the SSCHA calculation (2x2x2 in your case) you should obtain the same Hessian frequencies (regardless the centering is applied or not).

[tveskaeg]

Dear Raffaello,
thank you for your answer. I've retried the interpolation without ASR and centering (commands tensor3.Center() and tensor3.Apply_ASR()), but the frequencies vary only about < 0.1 cm^-1 compared to the interpolation including ASR and centering. So, the difference to the original Hessians is still much larger.

I'm still using CC.Spectral.get_static_correction_along_path for the interpolation.

All the best
Roman

[rafbianco]

Dear Roman,

I had a look at your script. I think I have found the problem. You are using as second order FCs the Hessian. This is not correct, because it is like doublecounting the odd correction. As second order FCs you must use the SSCHA dyn obtained at the end of the minimization.

However, this is not the first time that we see this misuse. We will be clearer in the tutorial about this point.
Thank you for your feedback.

Best,
Raffaello

[tveskaeg]

Great, thank you so much! This seems to work now. I didn't think too much about that since the comment in the first interpolation of the PbTe tutorial says "# For example, we can use PbTe.dyn#q, or PbTe.SSCHA.dyn#q, or PbTe.Hessian.dyn#q", so I just took the Hessians (although in the code lines you are correctly using *SSCHA.dyn#q).

I'll update this thread in the next days with the new figures for the higher q-grid interpolation!

Best,
Roman

[rafbianco]

Yes, I understand what is at the origin of the misunderstanding. Just to be clear: to initialize the 3rd order FCs you need the structural informations, that you can retreive from any dyn (SSCHA, Hessian, ...). But when you need the full information encoded in the dynamical matrices (frequencies and polarization vectors), then you must use the appropriate ones.

Best,
Raffaello

[tveskaeg]

Alright, makes sense.
I've checked the interpolation more carefully, and yes, the wrong initialization was the problem. I could reproduce the original Hessians and dispersions with the correct (final SSCHA) dyns, the method CC.Spectral.get_static_correction_along_path, and the QE postproc tool q2qstar.x.

In the same way, the interpolation to higher q-grids seems to work as well (Gamma, H, and N are the q-points shared explicitly by all three grids):

The black curves show the "check" (2x2x2 q/k-grid), red and green show the interpolation from my 2x2x2 q-grid to a 4x4x4 and 6x6x6 respectively with a 12x12x12 k-grid (where the Hessian frequencies are converged within < 1e-3 cm^-1).

Thanks a lot for the nice support, @rafbianco!

EDIT: I checked the dispersion "by hand" and interpolated directly on Gamma-H with CC.Spectral.get_static_correction_along_path (blue asterisks in the figure): this follows pretty much the original dispersion.
The problem is somewhere in (my application of) q2qstar.x; I found some wrong frequencies in there.

[anyi-mxh]

Dear Roman

I saw your discussion on the SSCHA forum regarding the interpolation of the Hessian matrix, and it seems that you’ve made great progress, successfully interpolating the Hessian matrix to a higher precision grid. I am currently facing a similar issue, and I was wondering if you could share the script you successfully implemented?

Best,
mxh

[tveskaeg]

Hey there,

for practical reasons in my workflows, the dyn.Interpolate routine turned out to be the most comfortable for preparing the files for EPW. We discussed this in the thread above only for the auxiliary matrices (dyn_end_population*), but from a technical/procedural perspective, you can also use it for the hessian matrices (replacing dyn_end_population* by dyn_hessian*). However, as discussed above, this is not the most precise way, and is only valid if you carefully converged the supercell size of your SSCHA calculations. I could do this with the help of force-fields (ML potentials).

dyn.Interpolate won't give you "new" or more precise information, it's just performing a Fourier interpolation based on the input Hessian matrices; i.e. the resulting dispersion (or force constants) obtained from the original Hessian matrices and the one obtained from the interpolated matrices are exactly the same. In principle, you could also do the interpolation in QE/matdyn if you provide the full list of q points in the grid (not just the irreducible points) and paste them together afterwards. Therefore, it's only advisable to use this approach if you are sure that your supercell size is well converged (which can be very hard with pure DFT calculations).

However, for a first test it's fine, if your harmonic phonons on a 222 q grid are not too different from your target grid (was that 555?). You can find the script in this answer above: https://github.com/orgs/SSCHAcode/discussions/49#discussioncomment-853580
You can choose to use it with or without harmonic support matrices, but always check the resulting dispersions and see if they look reasonable. If you don't use support matrices, make sure that the q-point numbering matches the numbering of the dvscf files (rename them by hand if necessary).

As Đorđe (Djordje) pointed out in the other thread, check if you are relaxing the structure within SSCHA (minimizer.minim_struct=true, if there are any free Wyckoff parameters) and re-calculate the scf and dvscf files for that specific structure before using them in EPW.

Unfortunately, I don't have a working script using the accurate routine CC.Spectral.get_static_correction_interpolated(dyn, tensor3, T, new_supercell, k_grid); and I don't if it is working in the current SSCHA version, as it had a bug in the 2021 version. Maybe the developers can comment on this?

Kind regards,
Roman