Push get_free_energy_hessian_dev() and environment

[ManexAlkorta]

Function under development. Employs a wave-vector classification of polarization-vector to write the expression of the free energy hessian matrix into smaller pieces. Those pieces have wave-vector q indices, where only the symmetry independent elements are computed (same for force-constants). It significantly outperforms previous function in required RAM and time in dense q-grids.

This pull request comes with a parallel pull request in cellconstructor.

Paralelization works with the old installation (old_setup.py).

[ManexAlkorta]

Now openmp_dep flags are added to the meson.build, and fortran parallelization works fine

[mesonepigreco]

Hi @ManexAlkorta ,
I think this implementation should have a test inside the testsuite that verifies that the old get_free_energy_hessian and the new get_free_energy_hessian_dev are consistent and give the same results, both with d3 only and also with d4.
We can even use a small 3x3x3 supercell of a gold tutorial (1 atom), which runs in CI only the d3 (but both offline, maybe tagging the test as a release only test) with an unconverged number of configurations (10?) so that it is very quick to run.

I think it is important to have a regression test that enforces correctness. I am suggesting a 3x3x3 supercell, as it is the smallest one with non-TRI q points (those gave me a lot of issues with the q-space Hessian in the TD-SCHA algorithm).

[mesonepigreco]

I see the addition here of OpenMP, but where is the OpenMP parallelization inside the Fortran subroutines? Maybe something is not clear to me.

[ManexAlkorta]

This is an issue I found with the new meson installation (old_setup.py works). It's not something about the new subroutines that I added, but was also failing working in parallel the old fortran subroutines (for example get_v3.f90). I use OpenMP is lots of subroutines inside module_hess.f90 (for example line 160 --> !$omp parallel private (i, nat1, nat2, nat3, alpha, beta, gamma, v3, indep_3fc, tmp_ref_3fc)).

[mesonepigreco]

I see that there are a lot of new variable defined in init().
This could be a problem as some of them are very heavy:

self.new_pol = np.zeros( (Nsc, Nsc * 3, 3), dtype = np.double)
self.ur = np.zeros( (self.N, Nsc * 3)) 
self.upsilon = np.zeros( (Nsc*3, Nsc * 3))

In particular, new_pol and upsilon are a materialization of a supercell dyn, since 1.4 we have removed all arrays Nsc x Nsc , as these are extremely slow with a big supercell (and memory-heavy).

Moreover, one of the problem we are facing on big ensembles is the storage of configurations and forces, and having an extra ur of dimension N_configs x Nsc*3 put a real burden on some SSCHA calculation (I have sscha calculations with big supercells of thousands of atoms, and tens of thousands of configurations where I cannot parallelize because multiple copies of the ensemble cannot be stored in memory, thus adding an extra one is heavy).

However, I see that these variables are only used in the new function get_free_energy_hessian_dev, so why not define them directly in that function? So a standard SSCHA minimization does not have to carry them in memory.

[ManexAlkorta]

Sure, I think this is the good way to go. Now all the new attribute initialization that I added are removed, and initialized inside get_free_energy_hessian_dev().

[ManexAlkorta]

I added a test of the hessian calculation a 4 atom SiC in a 3x3x1 q-grid (36 atoms). I would skip 1 atom systems for testing this function, as the most critical subtleties lay in the mapping and permutation of atoms/wave-vectors. This test takes around 10 minutes in my laptop, so i labeled it as a release test (not to be run at any minor change of the code).

[mesonepigreco]

This could be large depending on the system. especially if one must comparenit with the effect of the fourth order that could be smaller than the error, you have like a benchmark figure showing rhe sscha frequency vs hessian v4 old and hessian v4 new, to quantify the error relatively to the actual v4 induced shift?

[ManexAlkorta]

I performed these checks, and even though the v4 correction is not huge, it's enough to distinguish it from any error coming from the new method.

Zoom:

However it's important to note that this new method is not exact. It consinsts in a self-consistent convergence of $\Theta_{abcd}$ of Bianco et al (2017). The threshold for the scf loop to stop is a parameter of get_free_energy_hessian_dev(). However, I'll change the atol to atol=5e-5 (the actual max error is around 1e-5).

[mesonepigreco]

Ok, very good. The Theta_{abcd} has converged as a geometric series (I think Unai tried to implement something like that in the past)? One must be careful with geometric series, as I remember from Unai's work that in some systems the geometric series is divergent (usually extremely anharmonic and low-symmetry ones, like ice and high-pressure molecular hydrogen).

I see now there is a conflict with the main branch (probably because I just merged another PR). I have fixed it; wait for the test suite, and then we can merge.

[mesonepigreco]

Hi @ManexAlkorta, all tests passed. Very good. When you want me to merge, just comment, and we will merge. I am planning a release 1.7, probably in September (beginning of autumn), which I think will include our two new Hessian methods as the default.

[ManexAlkorta]

Great! From my side I think it's interesting to merge it now, so people can play with the function before the official release 1.7.

[mesonepigreco]

Ready to merge