[Q-e-developers] Loops on real-space grid

[Paolo Giannozzi]

On Mon, May 29, 2017 at 2:34 PM, Oliviero Andreussi
<oliviero.andreussi at usi.ch> wrote:

> Having said this, I double checked inside PW/src and the same patch could be done
> for PW/src/compute_dip.f90 and for PW/src/add_efield.f90

add_efield.f90 seems fine to me, thanks to this check:

     ! ... do not include points outside the physical range

     IF ( i >= dfftp%nr1 .OR. j >= dfftp%nr2 .OR. k >= dfftp%nr3 ) CYCLE

compute_dip.f90 accesses the charge density in nonphysical points and
sums over them, but, as long as the charge density is properly set to
0 in nonphysical points, this should not lead to wrong results. I
strongly prefer that nonphysical points are never accessed, though.

By the way: this problem of FFT dimensions has surfaced more than once
in the past; it is almost forgotten since the demise of AIX machines
with ESSL, because the trick of increasing first dimension by one is
not implemented for FFTW and Intel FFT. Not sure whether such trick
may still be useful in the future or whether it is a legacy of a
distant past.

> while add_monofield already had the test in place [...]

>      ! ... do not include points outside the physical range
>
>      IF ( k >= dfftp%nr3 ) CYCLE

I don't think this is sufficient, though: one should check the index
in xy plane as well.

> I could easily patch PW/src/compute_dip.f90 (which is, btw, very very similar
> to Modules/compute_dipole.f90, is there really a need for a duplicate?)

I have tried, and given up, more than once to merge these routines:
they look the same but there are subtle differences.

> In particular, it seems to me that while the loops in Modules have this general preliminary setup
>  [...]
> the same loops in PW/src have this slightly different one
>  [...]
> I understand that the two should be equivalent, provided that all the components of the fft
> descriptor are initialised correctly with and without MPI, so shall we only use one of the two
> approaches? Which?

I have a slight preference for the latter because it is the same for
MPI and non-MPI case. I would like to see #ifdef __MPI confined to
what is really different in serial and in parallel execution

Paolo

-- 
Paolo Giannozzi, Dip. Scienze Matematiche Informatiche e Fisiche,
Univ. Udine, via delle Scienze 208, 33100 Udine, Italy
Phone +39-0432-558216, fax +39-0432-558222