Ensemble sampling structures is generated too far away

[chenyubi14]

Dear SSCHA developers,

I'm encountering an issue with the ensemble generation that seems to be affecting my current results.

I am working on a CDW structure similar to NbSe2. When I calculate the low-temperature phonon dispersion (5K, 20K, 50K) by SSCHA, I find that, surprisingly, the high-symmetry pristine structure exhibits a stable phonon dispersion. I have followed the SSCHA procedure correctly, with 500 configurations for relaxation and 5000 configurations for Hessian. I suspect that the problem might be the ensemble generation, as it is sampling the random structures too far away. Below is the hessian phonon dispersion for 20K of the high-symm pristine structure.

My energy landscape for this material is below. The high-symmetry pristine structure has a higher energy than the low-symmetry CDW structure. This means I should expect an unstable phonon dispersion at very low temperatures like 20K. The distance between the pristine and CDW structure is around 5 [atomic_mass_unit^0.5 A]. (For simplicity, I will denote this distance unit as [A] below.) (You may ignore the different legends because these are the machine-learning force fields to approximate DFT.)

However, when I generate the ensemble sampling at 5K and 20K, the random structures are always > 10 [A] away from the high-symm structure. It means the sampling does not include the 5 [A] closer regions, which contain lower-energy CDW states. For structures far away, it will make the energy landscape look like a parabola, a local minima around high-symm pristine structure. I think this is why it gives me stable phonon dispersions.

I have tried to look at the ensemble generation code cellconstructor/Phonons.py. I am struggling to understand the implementation of ExtractRandomStructures, so I am reaching out here to seek help. In addition, I also found that 5K and 20K emsemble structures have similar distances, both are around 13[A] away, which means lower temperature is not going to result in an unstable phonon.

I would greatly appreciate any insights or guidance you could provide to help resolve this issue.
Thanks,
Yubi

[mesonepigreco]

Hi Yubi,
At very low temperatures it is normal that your system display very large displacements if the free energy landscape is shallow. Indeed 5A is a bit too much. The reason is that, likely, you have a very small frequency in your dynamical matrices, which results in the large displacements generated by the code

The equation is for the mean squared displacement is

$$\sigma = \frac 1m \sqrt{\frac{2n +1}{2\omega}}$$

So if you have a frequency that is very small, the code generates very big displacements.

The resulting sscha phonons are always positive, what decides if phonons are imaginary is the frequencies of the free energy hessian, calculated with the get_free_energy_hessian function.
Are you checking the hessian frequencies to inquire the stability of rhe structure?

All the best
Lorenzo

[chenyubi14]

Hi Prof. Dangic,

Are you sure that the unstable phonon mode is inside of your SSCHA supercell?

Yes, I am sure the CDW phonon mode is inside the SSCHA supercell size. With the same cell size, this was how I generated the energy landscape plot.

Just a comment on the sizes of displacements at 5K and 20K. You would expect them to be similar for your system. 20 K is about 0.5 THz. Whatever supercell you are using for SSCHA, you probably do not have any phonon modes with frequencies less than 0.5 THz, meaning a temperature increase from 5K to 20K does not affect your SSCHA procedure much.

Thanks for the comment. That makes sense to me. Then similar displacements between 5K and 20K are expected, but I need to figure out why 20K generates ensemble structures that are 13 [A] away. Do you have references for the formula for generation? Specifically, ExtractRandomStructures of cellconstructor/Phonons.py

a_mu = 1 / np.sqrt( np.tanh(beta*ws / 2) *2* ws) * BOHR_TO_ANGSTROM
total_coords = np.einsum("ij, i, j, kj->ik", pol_vects, 1/np.sqrt(mass1), a_mu, rand)

Thanks a lot!
Yubi

[DjordjeDangic]

Hello Yubi,

The formula for the mean square displacement that Lorenzo mentioned and that is implemented in the code is pretty much in every textbook that has a phonon section.

[mesonepigreco]

Dear Yubi.
That is a typical textbook equation. For example, in https://en.wikipedia.org/wiki/Quantum_harmonic_oscillator
you find the following equation:

It is quite straightforward to derive and most likely covered in any quantum statistical physics book.

However, from your dispersion, I do not see why you should have such large displacements. How are you computing the displacements? 20 A are sufficient to destroy completely the structure: If you visualize with a software your 2D structures, do you find crazy data? If so, it could be a problem of the machine learning force field that goes in a interpolation regime (where you have no more training data from the DFT) and converges on a fake structure.

[chenyubi14]

Dear Lorenzo and Đorđe,

Thanks so much for the response. I get your point. Here $a_{mu} = \sqrt{\langle x^2 \rangle}$ is the mean squre displacement. $\langle x^2 \rangle$ is a function of $\langle n \rangle$ which is a function of $\tanh(\beta \omega/2)$ through cannonical ensemble of harmonic oscillator. That makes a lot sense to me.

However, from your dispersion, I do not see why you should have such large displacements. How are you computing the displacements?

I compute the distance $Q$ by mass-weighted coordinate difference. $Q^2 = \sum_i M_i |R_{1,i} − R_{2,i}|^2$, where $M_i$ is mass of atom $i$; $R_1, R_2$ refers to the coordinates of two structures. It's worth noting that $Q$ scales linearly with the size of the supercell, doubling if the cell size doubles.

20 A are sufficient to destroy completely the structure: If you visualize with a software your 2D structures, do you find crazy data?

I've visualized some structures from the ensemble, and they don't appear significantly different from the high-symmetry pristine structure. I attribute this to the relatively large size of my supercell (4 * 4 * 2 of the unit cell, containing 288 atoms). My previous energy landscape plot was also calculated with 288 atoms.

Upon double-checking with other example materials like NbSe2 and silicon, I've observed that their ensemble structures also exhibit displacements of around 10 [A] with similar 300-atom cell sizes under a ultra-low temperature T=10K.

Therefore, correct me if I am wrong, it seems that achieving a 5 [A] displacement for a 288-atom cell size, particularly for the CDW material I am interested in, might be unattainable. This suggests that I could not produce the gradual disappearance of imaginary frequencies with finite-temperature calculations.

If so, it could be a problem of the machine learning force field that goes in a interpolation regime (where you have no more training data from the DFT) and converges on a fake structure.

I checked the machine learning force field (MLFF) works well for random structures. I have used some SSCHA generated structures as testing data, and MLFF works fine with energy, force, and stress data.

Thanks a lot for your helpful insights!
Yubi

[mesonepigreco]

Hi,
Ok, I see. So your displacements are not at all big, and 10 A is rescaled for the mass (which is high), generating a small displacement. Then I believe that the displacements you are getting are due to quantum fluctuations, as even for n=0, you still get a nonzero displacement, which is the quantum zero-point motion.

If this is the case, then your structure is stabilized at T=0K by quantum effects and the CDW does not occur.