H3S vc-relax weird behavior - Free Energy not decreasing?

[adenchfi]

Hello,
The SSCHA is an amazing idea! Thank you all for developing it.

I tried running the H3S example on the Bebop cluster at ANL, where my main change was running the calculations at fixed pressure instead of fixed volume, and changing the # of kpoints. I seem to get strange results:

The free energy has that strange jump, and then the code crashes with a segfault without any other information. I suspected it could be because of the Kong Liu effective sample size parameter getting smaller, so I tried to adapt the code to use the Relax module used for the LaH10 example. Here's the modified code:
H3S_initphon.txt
SSCHArelax.txt

The H3S_initphon script seems to work perfectly fine. The SSCHA relax script has the following output behavior:

The free energy has this strange behavior, and those large 'wide' regions are actually huge oscillations between steps. The code runs for about 6 iterations of the relax process, and then crashes.
It segfaults with the message /var/spool/slurmd/job1999230/slurm_script: line 23: 1000 Segmentation fault (core dumped) python3 SSCHArelax.py > relaxoutput.out. It doesn't appear to be a memory issue as it segfaults at the same spot whether or not I use more nodes.

If I change the # of kpoints, it changes how many iterations the code goes for before segfaulting. For 2x2x2 kpoints, it ran for 4 iterations before segfaulting.

The only thing I can think of is that the code isn't fully supporting of fixed pressure (as opposed to fixed-volume) calculations? But I'm not sure.

Looking forward to hearing back,
Adam

[mesonepigreco]

Hi, the first case of minimization is actually very good. The jump you see in the free energy is normal and simply a consequence of the fact you generated a new ensemble of ionic configurations. In this way, the free energy (which contains a term computed averaging the energy of the extracted configurations) will change in the new ensemble, but, as you see, remains compatible with the previous value within the stochastic error.

Actually, in the first case, I think the code converged (both gradients reached zero within their stochastic error). It is quite surprising it ended with a seg. fault. can you attach the full output so we can try to spot what went wrong?
Did the code print the dynamical matrix of the second population (should be named dyn_pop2_X)?

In the second case, if you optimize at a fixed pressure, you are minimizing the Gibbs free energy. The code is printing the Helmoltz free energy, so you should probably add a PV term to the free energy to see that one. Since it is a variable cell relaxation, your volume changes abruptly after each ensemble to optimize the Gibbs free energy. This is probably giving the jumps in the free energy. In particular, it seems you are optimizing the cell at a bigger pressure than the starting one, this means your volume is getting smaller after each population (to optimize the target pressure). Indeed, the Helmoltz free energy increases a lot when the volume decreases.
However check the pressure printed in output by the code:

cat output | grep "STRESS" -A 5

It will print the stress tensor after each minimization. You can check if it is correctly converging to the pressure value you choose.

Also in the second case, if you look at both the gradient and the Kong-Liu (Neff / N0) ratio, the code seems to be converged after almost every iteration (apart from the first). You recognize this because the Neff value reaches 1/2 of the initial value (N0) only in the first population. In this case the code halts because the ensemble is no more good, and a new ensemble is generated. However, in the next populations, the code stops when Neff > N0 / 2, this means that it stops because the gradient converged to 0 within the stochastic error (the simulation converged). You can check this looking at the gradient plot. Therefore, also here the code probably stopped because it reached the correct value of the pressure.
Again, it is a bit surprising this seg. fault error, If you can attach the output of relaxoutput.out we can check what are the last things the code does before crashing.
Bests,
Lorenzo

[adenchfi]

Hello,

Thank you for the thorough response! I learned a few things from it. Yes, I am getting "/var/spool/slurmd/job1999230/slurm_script: line 23: 1000 Segmentation fault (core dumped) python3 SSCHArelax.py > relaxoutput.out" at the end, after AMASS: {} {} is printed.

Extracting the pressure via the cat STRESS command, it appears the pressure I have started at 219 GPa (I tried relaxing it to 210 GPa), but it increases during each SSCHA step, up to 420 GPa or so. Hm. That is strange. My input desired pressure is 210 GPa, and it stopped once the code hit ~ 420 GPa. Maybe there's something accidentally doubling the target pressure?

Attached is the output text file for my 2nd case, along with input files for the phonon calc (H3S_initphon.py; this ran fine) ahd the SSCHA relaxation (SSCHArelax.py).

zipoutput.zip

Perhaps me specifying both 'press: 2100' in the SCF file input, and "target_press=2100" in the relax.vc_relax call resulted in a doubling to 420 GPa? While I wait for a reply, I'll try running things without the "press: 2100" in the SCF file input.

[mesonepigreco]

The most important files are nomm_spectral_func_[smear].dat, which stores the spectral functions. Then the code extracts from the spectral function useful info like the frequencies and lifetimes in the Lorentzian approximation. The most relevant info is, however, inside the spectral function. This bug seems to take place in this last part of the code, so it does not affect the spectral function itself. So:

In particular, is "nomm_spectral_func_[smear].dat" the final data file I want?

Yes!

@rafbianco @marcocherubini
Have you seen this issue before? I remember something like that, but I thought we fixed it.
Can you have a look at it?

[marcocherubini]

Hi,
It seems like the phonon frequencies you get are far outside the energy range you are setting.
If I understood correctly, you said that you are using the SSCHA free energy hessian dynamical matrices as the input for get_diag_dynamic_correction_along_path, but you should use the SSCHA dynamical matrices (the ones you get after the minimization).
In fact, in the hessian matrices you are already including the static correction. Basically, you are applying twice the correction
(even if the second time in the dynamical approximation).
This is for sure a source of error, even if I don't think it is the reason of the error you are getting. The static correction reduces the frequencies, so, it should not push them away from the energy range. Moreover, I think it is strange that in a bigger grid the bug disappears, since in any case the [1,1,2] contains the [1,1,1] where you have a problem. If you think the phonon frequencies you get are reasonable, you should increase the energy range (e0,e1) and this should solve the problem.

P.S. For what concerns the memory issue, we had something similar, but now it is fixed (in the master).
So, it is strange you cannot use a bigger k-grid, the memory cost is the same as for the [1,1,1].

[adenchfi]

Thanks you both!

Yes!
Excellent. I plotted the spectral function and it seems fine, very interesting to look at.

It seems like the phonon frequencies you get are far outside the energy range you are setting.

Yep, I've fixed that now.

In fact, in the hessian matrices you are already including the static correction. Basically, you are applying twice the correction

Okay, good to know I should use the SSCHA minimized quantities and not the free energy hessian. I'll need to review the SSCHA papers again now that I've actually used the software :-)

P.S. For what concerns the memory issue, we had something similar, but now it is fixed (in the master).

Got it. To get it working on our cluster I'm using the docker version (as per above discussion); does the current docker link reflect these changes?

Unrelated but, what format is the q-path file expecting? It doesn't look like crystal coordinates. I thought it was 2pi/alat, but I'm not sure if the XGX_path file in the PbTe example is reflecting that (or my own ignorance is again shining).

[rafbianco]

Hi,
the expected q-path format is 2pi/Angstrom.

[mesonepigreco]

Got it. To get it working on our cluster I'm using the docker version (as per above discussion); does the current docker link reflect these changes?

Not yet, I will prepare a new docker soon once we release the next alpha version of the code together with an update of the pip python packages. I am waiting to collect more bugfixes together.

@rafbianco

Okay, good to know I should use the SSCHA minimized quantities and not the free energy hessian.

Maybe we should change the tutorials on the spectral function specifying this more clearly with an alert, bold text or something like that, it seems a dangerous issue; something people could misunderstand.