[Q-e-developers] Loops on real-space grid
Oliviero Andreussi
oliviero.andreussi at usi.ch
Mon May 29 14:34:08 CEST 2017
Dear Thomas,
The bug was in a subroutine which exploited the real-space loop assuming that all the points were meaningful (fd_gradient, which computes the gradient of a function in the real-space grid using finite differences). In general, even though real-space loops exist in other subroutines in Modules and in PW/src, the test of wether the grid-point has physical sense or not is not necessary, because the operations inside the loop are on variables which have been already zeroed in those points. I added the tests to the routines in Modules because I didn't like too much having this sort of implicit assumptions in the code, but if it turns out that adding the ifs statements is inefficient from the performance point-of-view, I think even just adding some comments in the code would be ok.
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, while add_monofield already had the test in place (probably someone else figured out the problem before me):
"
! ... do not include points outside the physical range
IF ( k >= dfftp%nr3 ) CYCLE
"
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?), and PW/src/add_efield.f90, but before proceeding I wanted to have the opinion of the other developers, also concerning the standardised way of setting up these loops.
In particular, it seems to me that while the loops in Modules have this general preliminary setup
#if defined (__MPI)
index0 = dfftp%nr1x*dfftp%nr2x*SUM(dfftp%npp(1:me_bgrp))
ir_end = MIN(nnr,dfftp%nr1x*dfftp%nr2x*dfftp%npp(me_bgrp+1))
#else
index0 = 0
ir_end = nnr
#endif
!
DO ir = 1, ir_end
the same loops in PW/src have this slightly different one
idx0 = dfftp%nr1x*dfftp%nr2x*dfftp%ipp(me_bgrp+1)
!
DO ir = 1, dfftp%nr1x*dfftp%nr2x*dfftp%npl
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?
Best,
Oliviero
--
Senior Postdoctoral Researcher
École Polytechnique Fédérale de Lausanne (EPFL) and
Università della Svizzera Italiana (USI) of Lugano
USI Campus, Via G. Buffi 17, 6904 Lugano, Switzerland
Emails: oliviero.andreussi @ epfl.ch -or- usi.ch
Tel: +41-(0)58-666-4810 / Skype: olivieroandreussi
Web: https://sites.google.com/site/olivieroandreussi<https://sites.google.com/site/olivieroandreussi%E2%80%8B>
________________________________
From: q-e-commits-bounces at qe-forge.org <q-e-commits-bounces at qe-forge.org> on behalf of Thomas Brumme <thomas.brumme at uni-leipzig.de>
Sent: Monday, May 29, 2017 1:31:51 PM
To: Quantum ESPRESSO svn commit messages
Subject: Re: [Q-e-commits] r13536 - trunk/espresso/Modules
Dear Oliviero,
I don't know if the same bug occurs in cases using tefield (or
monofield) but those subroutines also rely on loops over real space.
Could you please add the tests also to add_efield.f90 and
add_monofield.f90 in PW/scr/
Thanks
Thomas
On 05/29/17 12:29, oliviero wrote:
> Author: oliviero
> Date: 2017-05-29 12:29:28 +0200 (Mon, 29 May 2017)
> New Revision: 13536
>
> Modified:
> trunk/espresso/Modules/compute_dipole.f90
> trunk/espresso/Modules/fd_gradient.f90
> trunk/espresso/Modules/generate_function.f90
> trunk/espresso/Modules/qmmm.f90
> Log:
> fixed a bug in fd_gradient (only called by Environ) and added some checks in other routines with loops on real-space grid.
>
> The bug was for the not-so-common case of machines where nr1x != nr1, loops in real-space should skip unphysical points.
>
>
> Modified: trunk/espresso/Modules/compute_dipole.f90
> ===================================================================
> --- trunk/espresso/Modules/compute_dipole.f90 2017-05-27 07:18:38 UTC (rev 13535)
> +++ trunk/espresso/Modules/compute_dipole.f90 2017-05-29 10:29:28 UTC (rev 13536)
> @@ -56,15 +56,18 @@
> ir_end = nnr
> #endif
> !
> - DO ir = 1, ir_end
> + DO ir = 1, ir_end
> !
> ! ... three dimensional indexes
> !
> i = index0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -82,7 +85,7 @@
> !
> CALL cryst_to_cart( 1, r, at, 1 )
> !
> - rhoir = rho( ir, 1 )
> + rhoir = rho( ir, 1 )
> !
> IF ( nspin == 2 ) rhoir = rhoir + rho(ir,2)
> !
>
> Modified: trunk/espresso/Modules/fd_gradient.f90
> ===================================================================
> --- trunk/espresso/Modules/fd_gradient.f90 2017-05-27 07:18:38 UTC (rev 13535)
> +++ trunk/espresso/Modules/fd_gradient.f90 2017-05-29 10:29:28 UTC (rev 13536)
> @@ -19,7 +19,7 @@
> USE kinds, ONLY: DP
>
> IMPLICIT NONE
> -
> +
> CONTAINS
> !=----------------------------------------------------------------------=!
> SUBROUTINE calc_fd_gradient( nfdpoint, icfd, ncfd, nnr, f, grad )
> @@ -39,7 +39,7 @@
> INTEGER, INTENT(IN) :: nnr
> REAL( DP ), DIMENSION( nnr ), INTENT(IN) :: f
> REAL( DP ), DIMENSION( 3, nnr ), INTENT(OUT) :: grad
> -
> +
> INTEGER :: index0, i, ir, ir_end, ipol, in
> INTEGER :: ix(-nfdpoint:nfdpoint),iy(-nfdpoint:nfdpoint),iz(-nfdpoint:nfdpoint)
> INTEGER :: ixc, iyc, izc, ixp, ixm, iyp, iym, izp, izm
> @@ -65,41 +65,44 @@
> #endif
> !
> DO ir = 1, ir_end
> - !
> + !
> i = index0 + ir - 1
> iz(0) = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( iz(0) .GE. dfftp%nr3 ) CYCLE ! if nr3x > nr3 skip unphysical part of the grid
> i = i - (dfftp%nr1x*dfftp%nr2x)*iz(0)
> iy(0) = i / dfftp%nr1x
> + IF ( iy(0) .GE. dfftp%nr2 ) CYCLE ! if nr2x > nr2 skip unphysical part of the grid
> ix(0) = i - dfftp%nr1x*iy(0)
> + IF ( ix(0) .GE. dfftp%nr1 ) CYCLE ! if nr1x > nr1 skip unphysical part of the grid
> !
> - DO in = 1, nfdpoint
> + DO in = 1, nfdpoint
> ix(in) = ix(in-1) + 1
> - IF( ix(in) .GT. dfftp%nr1x-1 ) ix(in) = 0
> + IF( ix(in) .GT. dfftp%nr1-1 ) ix(in) = 0
> ix(-in) = ix(-in+1) - 1
> - IF( ix(-in) .LT. 0 ) ix(-in) = dfftp%nr1x-1
> + IF( ix(-in) .LT. 0 ) ix(-in) = dfftp%nr1-1
> iy(in) = iy(in-1) + 1
> - IF( iy(in) .GT. dfftp%nr2x-1 ) iy(in) = 0
> + IF( iy(in) .GT. dfftp%nr2-1 ) iy(in) = 0
> iy(-in) = iy(-in+1) - 1
> - IF( iy(-in) .LT. 0 ) iy(-in) = dfftp%nr2x-1
> + IF( iy(-in) .LT. 0 ) iy(-in) = dfftp%nr2-1
> iz(in) = iz(in-1) + 1
> - IF( iz(in) .GT. dfftp%nr3x-1 ) iz(in) = 0
> + IF( iz(in) .GT. dfftp%nr3-1 ) iz(in) = 0
> iz(-in) = iz(-in+1) - 1
> - IF( iz(-in) .LT. 0 ) iz(-in) = dfftp%nr3x-1
> + IF( iz(-in) .LT. 0 ) iz(-in) = dfftp%nr3-1
> ENDDO
> !
> DO in = -nfdpoint, nfdpoint
> - i = ix(in) + iy(0) * dfftp%nr1x + iz(0) * dfftp%nr1x * dfftp%nr2x + 1
> - gradtmp(1,i) = gradtmp(1,i) - icfd(in)*f(ir)*dfftp%nr1x
> - i = ix(0) + iy(in) * dfftp%nr1x + iz(0) * dfftp%nr1x * dfftp%nr2x + 1
> - gradtmp(2,i) = gradtmp(2,i) - icfd(in)*f(ir)*dfftp%nr2x
> - i = ix(0) + iy(0) * dfftp%nr1x + iz(in) * dfftp%nr1x * dfftp%nr2x + 1
> - gradtmp(3,i) = gradtmp(3,i) - icfd(in)*f(ir)*dfftp%nr3x
> + i = ix(in) + iy(0) * dfftp%nr1x + iz(0) * dfftp%nr1x * dfftp%nr2x + 1
> + gradtmp(1,i) = gradtmp(1,i) - icfd(in)*f(ir)*dfftp%nr1
> + i = ix(0) + iy(in) * dfftp%nr1x + iz(0) * dfftp%nr1x * dfftp%nr2x + 1
> + gradtmp(2,i) = gradtmp(2,i) - icfd(in)*f(ir)*dfftp%nr2
> + i = ix(0) + iy(0) * dfftp%nr1x + iz(in) * dfftp%nr1x * dfftp%nr2x + 1
> + gradtmp(3,i) = gradtmp(3,i) - icfd(in)*f(ir)*dfftp%nr3
> ENDDO
> !
> ENDDO
> !
> #if defined (__MPI)
> - DO ipol = 1, 3
> + DO ipol = 1, 3
> CALL mp_sum( gradtmp(ipol,:), intra_bgrp_comm )
> CALL scatter_grid ( dfftp, gradtmp(ipol,:), grad(ipol,:) )
> ENDDO
> @@ -135,7 +138,7 @@
> !
> SELECT CASE ( ifdtype )
> !
> - CASE ( 1 )
> + CASE ( 1 )
> ! (2N+1)-point Central Differences
> IF ( nfdpoint .EQ. 1 ) THEN
> ncfd = 2
> @@ -153,12 +156,12 @@
> ncfd = 840
> icfd( 4 ) = -3
> icfd( 3 ) = 32
> - icfd( 2 ) =-168
> + icfd( 2 ) =-168
> icfd( 1 ) = 672
> ELSE
> WRITE(*,*)'ERROR: wrong number of points',nfdpoint,&
> &' for finite difference type ',ifdtype
> - STOP
> + STOP
> ENDIF
> !
> CASE ( 2 )
> @@ -166,12 +169,12 @@
> IF ( nfdpoint .GE. 2 ) THEN
> ncfd = (nfdpoint)*(nfdpoint+1)*(2*nfdpoint+1)/3
> DO in = 1,nfdpoint
> - icfd( in ) = in
> + icfd( in ) = in
> ENDDO
> ELSE
> WRITE(*,*)'ERROR: wrong number of points',nfdpoint,&
> &' for finite difference type ',ifdtype
> - STOP
> + STOP
> END IF
> !
> CASE ( 3 )
> @@ -185,19 +188,19 @@
> ncfd = 1188
> icfd( 4 ) = -86
> icfd( 3 ) = 142
> - icfd( 2 ) = 193
> + icfd( 2 ) = 193
> icfd( 1 ) = 126
> ELSE IF ( nfdpoint .EQ. 5 ) THEN
> ncfd = 5148
> icfd( 5 ) =-300
> icfd( 4 ) = 294
> icfd( 3 ) = 532
> - icfd( 2 ) = 503
> + icfd( 2 ) = 503
> icfd( 1 ) = 296
> ELSE
> WRITE(*,*)'ERROR: wrong number of points',nfdpoint,&
> &' for finite difference type ',ifdtype
> - STOP
> + STOP
> ENDIF
> !
> CASE ( 4 )
> @@ -215,19 +218,19 @@
> ncfd = 128
> icfd( 4 ) = 1
> icfd( 3 ) = 6
> - icfd( 2 ) = 14
> + icfd( 2 ) = 14
> icfd( 1 ) = 14
> ELSE IF ( nfdpoint .EQ. 5 ) THEN
> ncfd = 512
> icfd( 5 ) = 1
> icfd( 4 ) = 8
> icfd( 3 ) = 27
> - icfd( 2 ) = 48
> + icfd( 2 ) = 48
> icfd( 1 ) = 42
> ELSE
> WRITE(*,*)'ERROR: wrong number of points',nfdpoint,&
> &' for finite difference type ',ifdtype
> - STOP
> + STOP
> ENDIF
> !
> CASE ( 5 )
> @@ -241,19 +244,19 @@
> ncfd = 96
> icfd( 4 ) = -2
> icfd( 3 ) = -1
> - icfd( 2 ) = 16
> + icfd( 2 ) = 16
> icfd( 1 ) = 27
> ELSE IF ( nfdpoint .EQ. 5 ) THEN
> ncfd = 1536
> icfd( 5 ) = -11
> icfd( 4 ) = -32
> icfd( 3 ) = 39
> - icfd( 2 ) = 256
> + icfd( 2 ) = 256
> icfd( 1 ) = 322
> ELSE
> WRITE(*,*)'ERROR: wrong number of points',nfdpoint,&
> &' for finite difference type ',ifdtype
> - STOP
> + STOP
> ENDIF
> !
> CASE DEFAULT
>
> Modified: trunk/espresso/Modules/generate_function.f90
> ===================================================================
> --- trunk/espresso/Modules/generate_function.f90 2017-05-27 07:18:38 UTC (rev 13535)
> +++ trunk/espresso/Modules/generate_function.f90 2017-05-29 10:29:28 UTC (rev 13536)
> @@ -66,11 +66,16 @@
> !
> i = idx0 + ir - 1
> idx = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( idx .GE. dfftp%nr3 ) CYCLE
> IF ( axis .LT. 3 ) THEN
> i = i - (dfftp%nr1x*dfftp%nr2x)*idx
> idx = i / dfftp%nr1x
> + IF ( idx .GE. dfftp%nr2 ) CYCLE
> END IF
> - IF ( axis .EQ. 1 ) idx = i - dfftp%nr1x*idx
> + IF ( axis .EQ. 1 ) THEN
> + idx = i - dfftp%nr1x*idx
> + IF ( idx .GE. dfftp%nr1 ) CYCLE
> + END IF
> !
> idx = idx + 1 + shift
> !
> @@ -170,9 +175,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -283,9 +291,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -381,10 +392,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> - r = 0.D0
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -475,9 +488,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -572,9 +588,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -694,9 +713,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -789,9 +811,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
> @@ -860,9 +885,12 @@
> !
> i = idx0 + ir - 1
> k = i / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> i = i - (dfftp%nr1x*dfftp%nr2x)*k
> j = i / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> i = i - dfftp%nr1x*j
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> DO ip = 1, 3
> r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &
>
> Modified: trunk/espresso/Modules/qmmm.f90
> ===================================================================
> --- trunk/espresso/Modules/qmmm.f90 2017-05-27 07:18:38 UTC (rev 13535)
> +++ trunk/espresso/Modules/qmmm.f90 2017-05-29 10:29:28 UTC (rev 13536)
> @@ -487,10 +487,13 @@
> DO ir = 1, dfftp%nnr
> index = index0 + ir - 1
> k = index / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> index = index - (dfftp%nr1x*dfftp%nr2x)*k
> j = index / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> index = index - dfftp%nr1x*j
> i = index
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> s(1) = DBLE(i)/DBLE(dfftp%nr1)
> s(2) = DBLE(j)/DBLE(dfftp%nr2)
> @@ -633,10 +636,13 @@
> !
> index = index0 + ir - 1
> k = index / (dfftp%nr1x*dfftp%nr2x)
> + IF ( k .GE. dfftp%nr3 ) CYCLE
> index = index - (dfftp%nr1x*dfftp%nr2x)*k
> j = index / dfftp%nr1x
> + IF ( j .GE. dfftp%nr2 ) CYCLE
> index = index - dfftp%nr1x*j
> i = index
> + IF ( i .GE. dfftp%nr1 ) CYCLE
> !
> s(1) = DBLE(i)/DBLE(dfftp%nr1)
> s(2) = DBLE(j)/DBLE(dfftp%nr2)
>
> _______________________________________________
> Q-e-commits mailing list
> Q-e-commits at qe-forge.org
> http://qe-forge.org/mailman/listinfo/q-e-commits
>
--
Dr. rer. nat. Thomas Brumme
Wilhelm-Ostwald-Institute for Physical and Theoretical Chemistry
Leipzig University
Phillipp-Rosenthal-Strasse 31
04103 Leipzig
Tel: +49 (0)341 97 36456
email: thomas.brumme at uni-leipzig.de
_______________________________________________
Q-e-commits mailing list
Q-e-commits at qe-forge.org
http://qe-forge.org/mailman/listinfo/q-e-commits
More information about the developers
mailing list