Evaluate & Save System Energy from Force Constants

[ejmeitz]

Is it possible to enable a clean way to get the system energy evaluated within TDEP inside the extract_forceconstants command? I think this would be useful for assessing convergence of the force constants as the energy is what was used in the TDEP 3rd order paper to demonstrate convergence. Right now the energy is printed to the console or into the log file but that is not the easiest thing to parse and nearly impossible if you used MPI. Maybe I'm missing an alternate way to get this from TDEP if so please let me know! I do have my own code for this but it would be nice if to check convergence all I needed was TDEP and some basic python scripts. Thanks!

[flokno]

Do you mean the final model energy (for all samples) or the energy per sample?

[ejmeitz]

The energy for all of the samples used to calculate the force constants would be nice so you can calculate a MSE.

What is calculated internally? Based on the output below I assume you calculate the energy for all the samples already. The RMSE number is basically what I want as an output so I can make a plot of RMSE vs # of samples used to create the IFC and also RMSE vs r_cut.

ENERGIES (meV/atom):
rms std std(res)
input: 92.825302 5.284873 -
second order: 93.585004 5.502899 0.804262
third order: 0.921307 0.537397 0.587235

[flokno]

Yes, it is calculated. With --verbose you even get an output for some of the samples now. Using this flag or adding one to dump the energies to, say outfile.energies would be straightforward

[flokno]

I had one more look: The output from --verbose is broken for all-electron energies (like from FHI-aims), which is already fixed in my somewhat larger developer branch #16 . Maybe dumping all energies could be included there, since it already does the same thing for pressure (write pressure and stress per time step to file)

[ejmeitz]

Using this flag or adding one to dump the energies to, say outfile.energies would be straightforward

This would be perfect! If you do add it to one of the dev branches please let me know, I'd love to use it.

Does the current verbose flag give all the energies right now? I was not even aware this flag was a thing for extract_forceconstants. I'll have to check that out later.

[flokno]

No, it only prints a limited number of samples since the full list of samples can be in the thousands for the case of MD (which was the dominant way of sampling at the time)

[flokno]

Ok I added this in #38 . @ejmeitz please check that this is what you want and please double check if it is correct also. Note the unit which is given in meV/atom where no. of atoms is the no. of atoms in the supercell. That should be the same unit we use for the other output, also it is written to the header of the outfile.energies. That file should be as easy to read as np.loadtxt(...)

[ejmeitz]

Ok I added this in #38

I will test it this afternoon, thank you very much!

As a side note would you expect thousands of samples to be required for MD with an empirical potential (e.g SW silicon)? I was under the impression a couple hundred was usually enough, although like I said I have not checked my convergence.

[flokno]

As a side note would you expect thousands of samples to be required for MD with an empirical potential (e.g SW silicon)? I was under the impression a couple hundred was usually enough, although like I said I have not checked my convergence.

"enough" is a highly problem-specific word: https://github.com/tdep-developers/tdep-tutorials/tree/main/02_sampling/convergence

[ejmeitz]

@flokno this output looks perfect! Just to confirm E_pot is the value from my simulation (looks like it matches my dump.stat) and adding the other four columns (plus U0??) together gives the energy prediction from the taylor series using the force constants calculated by TDEP?

Edit: Nothing quite seems to make the energies add up to the MD energy I have. But the R^2 is 0.97

[flokno]

@flokno this output looks perfect! Just to confirm E_pot is the value from my simulation (looks like it matches my dump.stat) and adding the other four columns (plus U0??) together gives the energy prediction from the taylor series using the force constants calculated by TDEP?

Yes, this should be the case.

Edit: Nothing quite seems to make the energies add up to the MD energy I have. But the R^2 is 0.97

Did you include the baseline energy (U0)?

[flokno]

I checked and it seems work for me:

Things to note:

• The energy in infile.stat is in eV it needs to be scaled with the number of atoms in the supercell.
• If you use polar corrections, don't forget to add Epolar
• The baseline energy U0 is in eV/atom, not meV. I guess to be consistent it would be better to write the other energies in eV/atom as well, I have changed that in 9b93ac7

Does it work for you now?

[ejmeitz]

Yeah was a silly mistake I forgot to multiply the U0 by 1000. This looks good to me now.

While I have you here, what are the other numbers dumped in the energies.U0 file? My file has:
-0.415895270894E+01 -0.433074957112E+01 -0.432766705602E+01 -0.432766705602E+01
but based on the log info one of the last 2 is U0 I'm looking for.

[flokno]

Yes I see the problem. I made this more expliclity in 9b93ac7 as well. The output is now self-explanatory:

# Unit:      eV/atom
# no. atoms: 216
#            mean(Epot)                     mean(Epot - Epolar - E2)       mean(Epot - Epolar - E2 - E3)  mean(Epot - Epolar - E2 - E3 - E4)
            -0.273141372188E+05            -0.273141886234E+05            -0.273141886234E+05            -0.273141886234E+05
            ```