Some questions about the calculation of phonon free energy with dynaphopy

[Lweiweij]

Hello, I have encountered some problems when I calculate phonon free energy with dynaphopy, hope to get your help.

1. The harmonic phonon free energy calculated by dynaphopy and phonopy are very different, but both calculations use the same FORCE_SETS. I want to know whether dynaphopy calculates the harmonic phonon properties by calling phonopy? If so, why do I get the two results are different? This makes me can 't believe the result of the anharmonic phonon free energy calculated by dynaphopy.

2. I use the --thm_full tags to calculates the quasiparticle phonon free energy and the power spectrum phonon free energy. However, these results also are different, i also use the renormalized force constants to calculates anharmonic phonon free energy with phonopy, but this result is different from the previous two. Besides, By examining the mesh.yaml file generated by the phonopy calculation, I found that the phonon frequencies calculated using the renormalized force constant have eight zero frequencies in the Brillouin zone ( gamma point ). I want to know whether these situations means that my MD trajectory is not appropriate? My MD trajectory consist of 50000 time steps, relaxation of 500 ps before sampling, my time step is 0.001ps. Do I need to perform a longer simulation and sample a larger MD trajectory(100000 time steps)? By the way, in order to get a more accurate result, how much MD trajectory does dynaphopy usually need?

[abelcarreras]

Dynaphopy uses a mesh of 40 40 40 points to compute the thermal properties. Are you using the same mesh in phonopy. This may explain the differences in the thermal properties using the harmonic approach.

About the thermal properties results, dynaphopy calculates them using two different strategies.

1. quasiparticle: calculated from the normalized force constants.
2. Power spectrum: Calculated from the integration of the full power spectrum.

For method 1, the calculation of the thermal properties of quasiparticles requires a modification of the usual formula which is a little different from just the calculation of thermal properties using the renormalized force constants in phonopy. You can find the details here: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.064106

method 2 is intended to be used as a test because it requires a very good quality power spectrum and it will be quite sensitive to the quality of the power spectrum. In small highly symmetric unit cells it is usually quite close to quasiparticle, for large ones it might be more difficult to converge to the same value.

[Lweiweij]

Thank you for your reply!

1. I set the same mesh as the phonopy calculation in the input_file of Dynaphopy, So I think this may not be the reason for the different harmonic results.

2. According to your reply, should I use the quasiparticle free energy instead of power spectrum free energy for systems with many atoms(the unit cell i use have 400 atom)?

3. I also calculate the thermal properties using the MD trajectory of the output atomic velocity. This quasiparticle free energy is slightly different from the result obtained by using the MD trajectory of the output atom position. Is this small difference normal?

4. In order to get a more accurate result, how much MD trajectory does dynaphopy usually need?

[abelcarreras]

1. I will check then..
2. Yes, I recommend to use Quasiparticle free energy
3. You mean using velocities instead of positions as input? If this is the case, yeah, I would expect some small differences. dynaphopy computes the velocities from positions if the positions are provided so it would be slightly more accurate if you provide velocities directly.
4. This is very difficult to estimate since it is strongly dependent on your supercell (number of atoms, symmetry, etc..). You will need to check the convergence by yourself.

[Lweiweij]

Thank you for your quick reply!

1. I will wait for your advice on the first question at any time.
2. Yes, I used velocity and position as input respectively.

Among other Issues, I see what you said about the fact that a large change in the atomic configuration in the MD trajectory will result in not getting the correct results. But if my system is going to have a phase transition at a certain temperature (the position of some atoms will change very much), then in this case, can I still use these MD trajectories with a large change in structure as input to dynaphopy and get the correct result?
And because the system I study has an ordered-disordered phase transition, I can't get and use a structure of new phase (after the phase transition) to calculate the free energy.

[abelcarreras]

If your system goes though a phase transition then the phonon quasiparticle approximation will break at some point. All properties calculated from renormalized phonon frequencies (or force constants) rely in the fact that you have a small deviation respect to the harmonic model.

[abelcarreras]

I modified the code to make a little more clear the computation of the thermal properties. All this is currently available in the development branch if you want to test.

1. -thm tag. This computes the thermal properties from numerical integration of the DOS. This includes the quasiparticle corrections and it is the most reliable. For making sure this works one should check integration value. This indicates how good the numerical integration is (integral should be close to 1). This depends on setting a good frequency range in dynaphopy that covers all phonon frequencies. Harmonic here should be very close to phonopy thermal properties but contains some integration error. This should be used as reference to how well you can trust quasiparticle thermal properties.

For reference I included a new tag --thm_reference. This prints phonopy thermal properties calculations. Those are exactly what you get with phonopy and are printed in addition to the previous data. Quasiparticle properties here do not include the quasiparticle correction, just straightforward thermal properties of renormalized force constants.

2. -thm_full tag. This computed the thermal properties directly from the full spectrum (DoS) computed from the autocorrelation of the atomic velocities. This does not make any use of harmonic model/force constants. No quasiparticle approach, however requires a really huge supercell to have a straight good power spectrum sampling of all modes which in most cases is unfeasible.

Basically I recommend -thm using --thm_reference to compare with phonopy harmonic values to make sure that the integration is not very off. Also if needed it is possible to use --normalize_dos to rescale DOS to force the integration = 1.

If you try it, let me know how this works.

[Lweiweij]

I have submitted the task using the development version according to your suggestion(with "dynaphopy input_file xxx.lammpstrj -thm --thm_reference -ts 0.001 --temperature 300 --normalize_dos"), but my calculation may take some time (may be one day), and when it 's over, I will show you its results.

[Lweiweij]

I see that you mentioned “This depends on setting a good frequency range in dynaphopy that covers all phonon frequencies.”, do I need to add a tag similar to “-r” to set the range of the frequency?

[abelcarreras]

no, everything is controlled using -r

[Lweiweij]

The results still seem to have some problems.

[abelcarreras]

hmm, maybe the precision of the DOS is not good enough to compute the thermal properties in this case. Can your send me the harmonic force constants and unit cell so I can test by myself?
BTW, I can see that you used normalize_dos to rescale the DOS, what do you get in integration if you do not do that? is it far from 1?

[Lweiweij]

This is the result without using normalize_dos

Do you need MD trajectory? That file is very large.

[abelcarreras]

It seems that this small error in the integration produces a large difference in the free energy. I don't need the MD trajectory for now. I want to check which precision do I need in the DoS to compute the free energy reliably from it.

[Lweiweij]

OK, This Zip file contains the force constant and the unit cell.
test.zip

[abelcarreras]

What are the dimensions of your supercell? I am trying to generate the DOS from your FORCE_SETS using PHONOPY and I need this information?

[Lweiweij]

4x2x1

[abelcarreras]

I'm still testing it but meanwhile I just updated the code so now you can control the resolution of the DoS spectra used to compute the thermal properties. By default it uses the resolution of phonopy (as in previous versions) but setting --resolution flag now also controls the resolution of the DoS for thermal properties. This means that decreasing this value should improve the quality of the numerical integral used to compute the thermal properties and therefore the difference in the harmonic properties obtained by phonopy and dynaphopy should be smaller. If you want you can find this in the development branch.