EPH parameterisations

[Ruoyan0103]

Hi,

In the README file, I noticed a mention of an upcoming feature: “Create EPH parameterisations based on Nikita Medvedev's data.” I would like to ask whether it is possible to use Nikita Medvedev's data with the current model.

I found these Ce and kappa_e data from: here, and I will convert the units before I apply those data with the current model.

[artuuuro]

I should create an example tutorial for doing this.
Here are the steps:

1. Get the tables from Nikita Medvedevs site for the material you are interested in.
2. Tabulate them (and extrapolate if necessary) so that in the first column you have T in Kelvin and second column you have heat capacity in units of eV /Ang^3 / K and third column is heat conductivity in units of eV / ps / Ang / K
3. Drop the first column
4. add a header to the file with two columns in this form (the first value line corresponds to T = 0K

# comment line 1
# comment line 2
# comment line 3
N_values T_step
C_e[0] kappa_e[0]
C_e[1] kappa_e[1]
...
C_e[N_values-1] kappa_e[N_values-1]

I hope this helps. Write a follow-up if you need more help.

[Ruoyan0103]

Many thanks! I understand.

Actually, I encountered a problem with the data from that repository, where the electronic heat capacity C_e is very small.

When running the simulation with a very small C_e, around 1e-9, the system temperature rises uncontrollably. Only when I set C_e to 1e-6 does the simulation run normally.

Some paper mentioned that this is a numerical issue, and increase C_e manually 1. However, I checked the code, and it seems there is a routine, similar to the Courant–Friedrichs–Lewy (CFL) condition, that monitors numerical stability and adjusts the simulation time step.

If this is purely a numerical issue, then the CFL-like condition should handle very small C_e, right? I am wondering why the abnormal temperature increase still occurs.

[artuuuro]

I believe you are working on a semiconducting material where this is the Ce behaviour (based on that repository). The check in the code might not be sufficient for this system. There might be cases where Te becomes negative due to large fluctuations. You could try making the electronic system much larger than the ionic which should reduce fluctuations. But if you want to have temperature gradients in the system then this would not work. Increasing the Ce could help, but at the cost of physics. I am planning to improve the code for non-metallic systems in the near future (a year or so) and then I would look into this in more detail for now you will have to find a work-around.

…

On Wed, 24 Sept 2025, 09:04 Ruoyan0103, ***@***.***> wrote: *Ruoyan0103* left a comment (llnl/USER-EPH#12) <#12 (comment)> Many thanks! I understand. Actually, I encountered a problem with the data from that repository, where the electronic heat capacity C_e is very small. When running the simulation with a very small C_e, around 1e-9, the system temperature rises uncontrollably. Only when I set C_e to 1e-6 does the simulation run normally. Some paper mentioned that this is a numerical issue, and increase C_e manually 1 <https://doi.org/10.1016/j.nimb.2020.09.025>. However, I checked the code, and it seems there is a routine, similar to the Courant–Friedrichs–Lewy (CFL) condition, that monitors numerical stability and adjusts the simulation time step. However, even with this mechanism, the simulation becomes unstable with very small C_e, and I am wondering why this still happens. — Reply to this email directly, view it on GitHub <#12 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AANGQDPJYULW7PEFAMR2YMT3UI67FAVCNFSM6AAAAACGD6LNIKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTGMRWHA3TINBXGE> . You are receiving this because you commented.Message ID: ***@***.***>

[Ruoyan0103]

Many thanks.

You could try making the electronic system much larger than the ionic which
should reduce fluctuations. But if you want to have temperature gradients
in the system then this would not work.

I tried this method (where the electronic system is 20 times larger than the atomic system), but I also had to reduce the timestep to 0.01 fs; otherwise, the temperature still rises too high.

Increasing the Ce could help, but at the cost of physics.

Instead of increasing Ce, I also tried increasing the rho value from 1 to 1000, which resolved the current problem. However, I am not sure how this will affect the physics.

[correaa]

Increasing rho would be equivalent to increasing the electron-ion coupling (if we are talking about the same rho), that is physics too.

There is yet another option that doesn't affect the physics that is to decrease the time-step, which if you can afford it, is the best option. (I say afford because the simulations will be slower).
The PDE stability condition is not exact, is only a local approximation, the true value it may depends on boundary conditions and maybe "what the ions are doing". So there could be a lower multiplicative factor that needs to be used in this case.

If what you are trying to model a semiconductor, take a few things into account. First, semiconducators have more complicated electron populations than a metal, electrons and holes can be at different temperatures and sustain different chemical potentials.

Although this is physics that will be missing here, I think it is still worth simulating semiconductors because it might not be an important factor that affects the ion dynamics, which is what we are ultimately doing.

Also, take into account that Ce -> 0 not only problematic numerically, but also probably unphysical for a semiconductor, first because deformations of the perfect lattice can effectively lower (or close) the gap that is already small (increasing the Ce).
For example, Ce -> 0 could be an artifact of calculating Ce from static DFT energies for a perfect crystal as a function of (electronic) temperature.

So, another option is to have a "floor" for the Ce values at low temperature, to workaround this.
For example, DFT's Ce with ions that are not exactly at equilibrium positions will be generally higher.

Another "physical solution" is to start the system with temperatures not lower than the Debye temperature in the first place. The simulated system is still the "classical" limit of the electron-phonon coupling, and that includes the phonons by themselves, so it doesn't make sense to simulate the system below the Debye temperature anyway.

To summarize, these are some (complementary) strategies,

1. lower the time-step (of the multiplicative factor of the time-step that comes from the stability condition). (given your illustrative experiment with Ce values, it is possible that the prefactor is missed by factor of sqrt(1000) ~ 30 but can be lower really.
2. start the simulation at a higher temperature (as high as you can afford for the system you want to mode), where Ce is not so small, and/or no lower than the Debye temperature.
3. if all fails, come up with a floor value for Ce, if that floor value is physical, better yet (e.g. from fluctuations of the ionic positions, either from temperature or from Zero-point energy).

I hope it helps.