Fail to reproduce the phonon dispersion of CDW materials at low temperatures

[hellolori]

Dear developers,

I have recently been studying the use of SSCHA combined with Quantum Espresso to calculate the phonon spectrum of CDW materials at low temperatures. However, I have encountered a peculiar issue. For single-layer CDW materials such as NbSe2, NbS2, and VSe2, the phonon spectrum I calculated using SSCHA at zero temperature no longer exhibits imaginary frequencies, which is inconsistent with the results reported in the literature. I will use the calculation of single-layer NbSe2 as an example and hope that you can help me some advices.

The phonon spectrum of this system has been reported in the literature [Phys. Rev. Lett. 2020, 125: 106101] and [J. Phys. Condens. Matter 2021, 33: 363001]. I adopted the same computational parameters as those used in [Phys. Rev. Lett. 2020, 125: 106101] and was able to obtain a consistent harmonic phonon spectrum, but I could not reproduce the anharmonic phonon spectrum at 0 Kelvin.

The left panel in the figure below shows my calculation results (using a 6x6x1 supercell), while the right panel is Figure 6 from the supplementary material of [Phys. Rev. Lett. 2020, 125: 106101]. As can be seen, in the literature, imaginary frequencies still exist at 0 Kelvin (purple solid line).

Here are some intermediate data from my calculation process. Judging from the change in the Effective sample size with the minimization steps, it appears that the calculation results have converged.

In the SSCHA calculation process, I adopted a manual approach. After completing the calculation of one generation of a structural ensemble, if the results do not converge, I generate the next generation of ensemble configurations based on the "final_sscha_dyn_population${index_popu}${iq}" file obtained from the last ensemble of the current generation, until the calculation converges. Finally, I use the "final_sscha_dyn_population${index_popu}${iq}" file from the last generation to plot the phonon spectrum.

Below is the calculation process for completing one generation of structural ensemble calculations:
(1) Generate a new generation of structural ensemble based on the dynamical matrix file obtained from the previous generation ensemble calculation.

import cellconstructor as CC, cellconstructor.Phonons
import sscha, sscha.Ensemble
import numpy as np
# Fix the seed so that we all generate the same ensemble
#np.random.seed(0)
# Load the dynamical matrix
dyn = CC.Phonons.Phonons("final_sscha_dyn_population4_", nqirr=7)
#
#[ apply here the needed changes to dyn ]
dyn.Symmetrize()
dyn.ForcePositiveDefinite()
#
# Prepare the ensemble
temperature = 0.0 #  K
ensemble = sscha.Ensemble.Ensemble(dyn, temperature)
#
# Generate the ensemble
number_of_configurations = 200
ensemble.generate(number_of_configurations)
#
# Save the ensemble into a directory
save_directory = "data"
population_id = 5
ensemble.save(save_directory, population_id)

(2) Combine the structure files with the QE header file "espresso.pwi" to generate input files for self-consistent calculations for each structure. The content of the espresso.pwi file is as follows:

&control
   prefix = 'nbse2'
   calculation = 'scf'
   restart_mode = 'from_scratch'
    tstress = .true.
   tprnfor = .true.
   disk_io='none'
pseudo_dir='/home/hellolori/pp'
/
&system
   ibrav = 0
   nat = 108
  ntyp = 2
   ecutwfc = 35
   ecutrho=400
   occupations  =  'smearing' ,
               smearing  =  'mp' ,
               degauss  =  0.005 ,
/
&electrons
   electron_maxstep = 150
   conv_thr = 1.0d-8
   mixing_beta = 0.3d0
   mixing_mode = 'plain'
/
ATOMIC_SPECIES
  Nb   92.91  Nb.pbe-spn-rrkjus_psl.1.0.0.UPF
  Se   78.96  Se_ONCV_PBE-1.1.upf
K_POINTS automatic
5  5  1 0 0 0

(3) Perform DFT self-consistent calculations for all configurations to obtain the total energy, atomic forces, and stress.

(4) Based on the data obtained from the DFT calculations, use SSCHA to minimize the free energy. The corresponding input file is:

import sscha, sscha.Ensemble, sscha.SchaMinimizer
import sscha.Utilities
import cellconstructor as CC, cellconstructor.Phonons
#
POPULATION = 5
#
dyn = CC.Phonons.Phonons("final_sscha_dyn_population4_", 7)
dyn.Symmetrize()
dyn.ForcePositiveDefinite()
temperature = 0.0 #  K
ensemble = sscha.Ensemble.Ensemble(dyn, temperature)
ensemble.load("data", population = POPULATION, N = 200)
#
minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minim.meaningful_factor=1e-5
minim.set_minimization_step(0.01)
minim.init()
#
# Save the minimization details
ioinfo = sscha.Utilities.IOInfo()
ioinfo.SetupSaving("minim_{}".format(POPULATION))
#
minim.run(custom_function_post = ioinfo.CFP_SaveAll)
minim.finalize()
minim.dyn.save_qe("final_sscha_dyn_population{}_".format(POPULATION))

(5) If step 4 does not converge, return to step 1 and generate the next generation of structural ensemble. Repeat the above process until step 4 converges.

Best regards,
hellolori

[DjordjeDangic]

In that reference, they are plotting the Hessian of free energy that has an imaginary phonon frequency. You are plotting, if I am not mistaken, the auxiliary phonons. You should compare your results with the dashed purple line from the reference paper and I would say you are spot on.

[hellolori]

Dear DjordjeDangic,

Thank you for your response. I indeed notice that my calculation results closely match the purple dashed line in the literature (right panel in the figure). However, these results are not reasonable because the CDW transition temperature of monolayer NbSe2 is above 100 Kelvin. Therefore, at 0 Kelvin, the system's phonon spectrum should exhibit imaginary frequencies, which correspond to the purple solid line.

My next question is whether it is possible to use the SSCHA to calculate the results in the purple solid line, and if so, how to proceed.

I noticed that the purple dashed line is labeled with $\Phi^{(3)}=0$, while the purple solid line has no such label. If my understanding is correct, the calculation results for the purple solid line take into account corrections from the third-order and higher-order force constant matrices, while the dashed line does not. Could you please explain how to specify the consideration of third-order or fourth-order force constant matrix corrections in SSCHA?

Thank you in advance for your guidance.

[DjordjeDangic]

Exactly, the auxiliary phonon frequencies that you are plotting right now can not be imaginary (if they are it is an interpolation problem). So they can not tell you anything about the phase transition.

Free energy Hessian is the quantity that can tell you at what conditions phase transition happens. To calculate free energy Hessian follow this tutorial: https://sscha.eu/Tutorials/tutorial_03_secondorder_phase_transitions/ . The function that will produce these free energy Hessian force constants and that will give you the full line from the reference is this one:

dyn_hessian = ensemble.get_free_energy_hessian(include_v4 = False)

Keep in mind that converging results for free energy Hessian usually requires a much larger number of configurations compared to the usual minimization procedure!

We also have video tutorials from the SSCHA school that took place last year (https://sscha.eu/lectures/). I strongly suggest you watch these because they go in-depth about the issue that you are interested in (especially the one by Raffaelo on second-order phase transitions).

[hellolori]

Thank you very much for your suggestions. I will first study the "tutorial_03_secondorder_phase_transitions/" section. After that, I will attempt to reproduce the calculation results reported in the literature again.