merge with the latest develop

CPV/src/wave.f90

@@ -88,21 +88,26 @@ subroutine elec_fakekine_x( ekincm, ema0bg, emass, c0, cm, ngw, n, noff, delt )
     real(DP) :: ftmp
     integer  :: i
 
-    ALLOCATE( emainv( ngw ) )
-    emainv = 1.0d0 / ema0bg
-    ftmp = 1.0d0
-    if( gstart == 2 ) ftmp = 0.5d0
-
     ekincm=0.0d0
-    do i = noff, n + noff - 1
-      ekincm = ekincm + 2.0d0 * wave_speed2( c0(:,i), cm(:,i), emainv, ftmp )
-    end do
-    ekincm = ekincm * emass / ( delt * delt )
+
+    IF( ngw > 0 ) THEN
+
+       ALLOCATE( emainv( ngw ) )
+       emainv = 1.0d0 / ema0bg
+       ftmp = 1.0d0
+       if( gstart == 2 ) ftmp = 0.5d0
+
+       do i = noff, n + noff - 1
+         ekincm = ekincm + 2.0d0 * wave_speed2( c0(:,i), cm(:,i), emainv, ftmp )
+       end do
+       ekincm = ekincm * emass / ( delt * delt )
+       DEALLOCATE( emainv )
+
+    END IF
 
     CALL mp_sum( ekincm, intra_bgrp_comm )
     IF( nbgrp > 1 ) &
        CALL mp_sum( ekincm, inter_bgrp_comm )
-    DEALLOCATE( emainv )
 
     return
   end subroutine elec_fakekine_x

FFTXlib/fft_helper_subroutines.f90

@@ -21,7 +21,7 @@ SUBROUTINE tg_reduce_rho_1( rhos, tg_rho_nc, tg_rho, ispin, noncolin, domag, des
      REAL(DP), INTENT(INOUT)  :: tg_rho_nc(:,:)
      REAL(DP), INTENT(OUT) :: rhos(:,:)
 
-     INTEGER :: ierr, ioff, idx, ir3, ir, ipol, ioff_tg, nxyp, npol_
+     INTEGER :: ierr, ioff, ir3, ir, ipol, ioff_tg, nxyp, npol_
 !     write (*,*) ' enter tg_reduce_rho_1'
 
 #if defined(__MPI)
@@ -71,7 +71,7 @@ SUBROUTINE tg_reduce_rho_2( rhos, tmp_rhos, ispin, desc )
      REAL(DP), INTENT(INOUT)  :: tmp_rhos(:)
      REAL(DP), INTENT(OUT) :: rhos(:,:)
 
-     INTEGER :: ierr, ioff, idx, ir3, nxyp, ioff_tg
+     INTEGER :: ierr, ioff, ir3, nxyp, ioff_tg
 !     write (*,*) ' enter tg_reduce_rho_2'
 
      IF ( desc%nproc2 > 1 ) THEN

FFTXlib/fft_types.f90

@@ -140,7 +140,7 @@ MODULE fft_types
 
   PUBLIC :: fft_type_descriptor, fft_type_init
   PUBLIC :: fft_type_allocate, fft_type_deallocate
-  PUBLIC :: fft_stick_index
+  PUBLIC :: fft_stick_index, fft_index_to_3d
 
 CONTAINS
 
@@ -994,6 +994,41 @@ PURE FUNCTION fft_stick_index( desc, i, j )
       fft_stick_index = desc%isind ( mc ) 
    END FUNCTION
 
+   !
+   SUBROUTINE fft_index_to_3d (ir, dfft, i,j,k, offrange)
+     !
+     !! returns indices i,j,k yielding the position of grid point ir
+     !! in the real-space FFT grid described by descriptor dfft:
+     !!    r(:,ir)= i*tau(:,1)/n1 + j*tau(:,2)/n2 + k*tau(:,3)/n3
+     !
+     IMPLICIT NONE
+     INTEGER, INTENT(IN) :: ir
+     !! point in the FFT real-space grid
+     TYPE(fft_type_descriptor), INTENT(IN) :: dfft
+     !! descriptor for the FFT grid
+     INTEGER, INTENT(OUT) :: i
+     !! (i,j,k) corresponding to grid point ir
+     INTEGER, INTENT(OUT) :: j
+     !! (i,j,k) corresponding to grid point ir
+     INTEGER, INTENT(OUT) :: k
+     !! (i,j,k) corresponding to grid point ir
+     LOGICAL, INTENT(OUT) :: offrange
+     !! true if computed i,j,k lie outside the physical range of values
+     !
+     i     = ir - 1
+     k     = i / (dfft%nr1x*dfft%my_nr2p)
+     i     = i - (dfft%nr1x*dfft%my_nr2p) * k
+     j     = i /  dfft%nr1x 
+     i     = i -  dfft%nr1x * j
+     j     = j + dfft%my_i0r2p
+     k     = k + dfft%my_i0r3p
+     !
+     offrange = (i < 0 .OR. i >= dfft%nr1 ) .OR. &
+          (j < 0 .OR. j >= dfft%nr2 ) .OR. &
+          (k < 0 .OR. k >= dfft%nr3 )
+     !
+   END SUBROUTINE fft_index_to_3d
+
 !=----------------------------------------------------------------------------=!
 END MODULE fft_types
 !=----------------------------------------------------------------------------=!

GWW/head/head.f90

@@ -104,7 +104,6 @@ PROGRAM head
      !
      !  If necessary the bands are recalculated
      !
-     !IF (setup_pw) CALL run_pwscf(do_band)
      IF (setup_pw) CALL run_nscf(do_band, iq)
      !
      !  Initialize the quantities which do not depend on

Modules/compute_dipole.f90

@@ -16,6 +16,7 @@ SUBROUTINE compute_dipole( nnr, rho, r0, dipole, quadrupole )
   USE kinds,            ONLY : DP
   USE cell_base,        ONLY : at, bg, alat, omega
   USE fft_base,         ONLY : dfftp
+  USE fft_types,        ONLY : fft_index_to_3d
   USE mp_bands,         ONLY : intra_bgrp_comm
   USE mp,               ONLY : mp_sum
   !
@@ -33,8 +34,9 @@ SUBROUTINE compute_dipole( nnr, rho, r0, dipole, quadrupole )
   ! ... Local variables
   !
   REAL(DP) :: r(3), rhoir
-  INTEGER  :: i, j, k, ip, ir, ir_end, idx, j0, k0
+  INTEGER  :: i, j, k, ip, ir, ir_end
   REAL(DP) :: inv_nr1, inv_nr2, inv_nr3
+  LOGICAL  :: offrange
   !
   ! ... Initialization
   !
@@ -46,28 +48,17 @@ SUBROUTINE compute_dipole( nnr, rho, r0, dipole, quadrupole )
   quadrupole(:) = 0.D0
   !
 #if defined (__MPI)
-  j0 = dfftp%my_i0r2p ; k0 = dfftp%my_i0r3p
   ir_end = MIN(nnr,dfftp%nr1x*dfftp%my_nr2p*dfftp%my_nr3p)
 #else
-  j0 = 0 ; k0 = 0
   ir_end = nnr
 #endif
   !
   DO ir = 1, ir_end
      !
      ! ... three dimensional indexes
      !
-     idx = ir -1
-     k   = idx / (dfftp%nr1x*dfftp%my_nr2p)
-     idx = idx - (dfftp%nr1x*dfftp%my_nr2p)*k
-     k   = k + k0
-     IF ( k .GE. dfftp%nr3 ) CYCLE
-     j   = idx / dfftp%nr1x
-     idx = idx - dfftp%nr1x * j
-     j   = j + j0
-     IF ( j .GE. dfftp%nr2 ) CYCLE
-     i   = idx
-     IF ( i .GE. dfftp%nr1 ) CYCLE
+     CALL fft_index_to_3d (ir, dfftp, i,j,k, offrange)
+     IF ( offrange ) CYCLE
      !
      DO ip = 1, 3
         r(ip) = DBLE( i )*inv_nr1*at(ip,1) + &

Modules/fd_gradient.f90

@@ -27,6 +27,7 @@ SUBROUTINE calc_fd_gradient( nfdpoint, icfd, ncfd, nnr, f, grad )
   USE kinds,         ONLY : DP
   USE cell_base,     ONLY : at, bg, alat
   USE fft_base,      ONLY : dfftp
+  USE fft_types,     ONLY : fft_index_to_3d
   USE scatter_mod,   ONLY : scatter_grid
   USE mp,            ONLY : mp_sum
   USE mp_bands,      ONLY : me_bgrp, intra_bgrp_comm
@@ -40,37 +41,27 @@ SUBROUTINE calc_fd_gradient( nfdpoint, icfd, ncfd, nnr, f, grad )
   REAL( DP ), DIMENSION( nnr ), INTENT(IN) :: f
   REAL( DP ), DIMENSION( 3, nnr ), INTENT(OUT) :: grad
 
-  INTEGER :: idx, j0, k0, i, ir, ir_end, ipol, in
+  INTEGER :: 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
   REAL( DP ), DIMENSION( :, : ), ALLOCATABLE :: gradtmp, gradaux
+  LOGICAL  :: offrange
   !
   grad = 0.D0
   !
   ALLOCATE( gradtmp( dfftp%nr1x*dfftp%nr2x*dfftp%nr3x, 3 ) )
   gradtmp = 0.D0
   !
 #if defined (__MPI)
-  j0 = dfftp%my_i0r2p ; k0 = dfftp%my_i0r3p
   ir_end = MIN(nnr,dfftp%nr1x*dfftp%my_nr2p*dfftp%my_nr3p)
 #else
-  j0 = 0 ; k0 = 0
   ir_end = nnr
 #endif
   !
   DO ir = 1, ir_end
     !
-    idx   = ir - 1
-    iz(0) = idx / (dfftp%nr1x*dfftp%my_nr2p)
-    idx   = idx - (dfftp%nr1x*dfftp%my_nr2p)*iz(0)
-    iz(0) = iz(0) + k0
-    IF ( iz(0) .GE. dfftp%nr3 ) CYCLE ! if nr3x > nr3 skip unphysical part of the grid
-    iy(0) = idx / dfftp%nr1x
-    idx   = idx - dfftp%nr1x*iy(0)
-    iy(0) = iy(0) + j0
-    IF ( iy(0) .GE. dfftp%nr2 ) CYCLE ! if nr2x > nr2 skip unphysical part of the grid
-    ix(0) = idx 
-    IF ( ix(0) .GE. dfftp%nr1 ) CYCLE ! if nr1x > nr1 skip unphysical part of the grid
+    CALL fft_index_to_3d (ir, dfftp, ix(0), iy(0), iz(0), offrange)
+    IF ( offrange ) CYCLE
     !
     DO in = 1, nfdpoint
       ix(in) = ix(in-1) + 1