issues with ph1d in energy.F90 and scfcv.F90 [SOLVED]

[sipr]

For some builds and/or machines, problems may occur in energy.F90 when ph1d is copied to gs_hamk%ph1d via
gs_hamk%ph1d(:,:) = ph1d(:,:)

Namely, the array ph1d is not defined for all indices (it is set in getphd called by scfcv), because its dimentions are determined by the largest of numbers ngfft(1),ngfft(2),ngff(3) while the "real" array is defined only up to a boundary determined by the actual values of the ngfft(1),ngfft(2),ngff(3) numbers.

As a result, when copying ph1d(:,:) to gs_hamk%ph1d(:,:) in energy.F90 , the code will touch also "undefined" values which have been initiated randomly and this may cause "floating point" errors, which are de facto irrelevant but which still make the code to crash.

A quick fix I am using is to defined the array ph1d in the scfcv.F90 subroutine for all indices in the beginning of this subroutine:
ph1d(:,:) = 1.0_dp ! Otherwise issues in energy due to undefined values

A more "grown-up" way of dealing with this might be either amending the dimension of the ph1d array or restricting the indices for which ph1d is copied to gs_hamk%ph1d in energy.F90. However, that would be too advanced for my degree of familiarity with the package.

I encountered this issue with the 6.12.1 version.
I understand that this error does not occur always and everywhere but when it occurs, it hurts...

Ondrej

[torrent]

Dear Ondrej,

Thanks for the bug report;
after several months, I have included a modification in Abinit (to be published in Abinit v7).
I have initialized ph1d directly in the getph subroutine (dedicated to ph1d computation).

Marc Torrent