Fortran interface: xcfun_eval_vec uses wrong pitch values

[AndHessel]

I think that the current pitch values

d_pitch = int(size(density) - 1, kind=c_int)
r_pitch = int(size(res) - 1, kind=c_int)

computed in xcfun_eval_vec are incorrect. For example, size(density) is the total length of the density array, not the size of the number of density values per grid point. In my opinion it should read

d_pitch = int(size(density(:,1)), kind=c_int)
r_pitch = int(size(res(:,1)), kind=c_int)

which appears to give consistent results when I use a batch of grid points that all have the same density values. In fact, with the original settings for d/r_pitch my program crashes when xcfun_eval_vec is called (either internally or sometimes on exit due to a problem with the release of memory).

[robertodr]

Thanks for reporting! Yes, it appears you're right. The xcfun_eval_vec function is tested only in the Python bindings and it could be that the same dimensioning mistake lurks there too. @chjacob-tubs @aspgomes could you have a look too? Would you mind opening a pull request with the proposed changes? It would be great if you could modify the example/Fortran_host/example.f90 to also test xcfun_eval_vec.

[AndHessel]

Ok, I can take a look (note that I have more experience with GitLab, but it should be not that much different to do merge requeste with GitHub, I guess). Note that I have not created Python bindings, so have not checked whether the bug appears in there, too.

Regarding an example for example.f90, I have written a small program to compare the results (pbex, order=2) against libxc. I can maybe extract the xcfun part of it and include it as another check in example.f90.

[chjacob-tubs]

I think in the Python interface this is handled correctly, see python/xcfun_export.cpp, lines 151-156:

            xcfun::xcfun_eval_vec(fun,
                                  nr_points,
                                  density.data(),
                                  density.shape(dens_ndim - 1),
                                  output.mutable_data(),
                                  output.shape(output_ndim - 1));

(I didn't write this code and I am not really familiar with PyBind11, but this looks good to me - and it is working correctly also in cases that there is more than one density value per point, which we test for the Python interface)

[robertodr]

Closed by #149. Thanks again for reporting the issue and providing a fix!