Avoid anything that any compiler may interpret as an attempt to access arrays

with zero dimension - GPU are especially hostile to zero-dimensional arrays

PW/src/allocate_nlpot.f90

@@ -52,7 +52,6 @@ SUBROUTINE allocate_nlpot
   ALLOCATE( nhtolm(nhm,nsp) )
   ALLOCATE( nhtoj(nhm,nsp)  )
   ALLOCATE( ijtoh(nhm,nhm,nsp) )
-  ALLOCATE( indv_ijkb0(nat)    )
   ALLOCATE( deeq(nhm,nhm,nat,nspin) )
   IF ( noncolin ) THEN
      ALLOCATE( deeq_nc(nhm,nhm,nat,nspin) )
@@ -66,6 +65,9 @@ SUBROUTINE allocate_nlpot
   ELSE
     ALLOCATE( dvan(nhm,nhm,nsp) )
   ENDIF
+  ALLOCATE (becsum( nhm * (nhm + 1)/2, nat, nspin))
+  IF (tqr) ALLOCATE (ebecsum( nhm * (nhm + 1)/2, nat, nspin))
+  ALLOCATE( indv_ijkb0(nat)    )
   ! GIPAW needs a slighly larger q-space interpolation for quantities calculated
   ! at k+q_gipaw, and I'm using the spline_ps=.true. flag to signal that
   IF ( spline_ps .AND. cell_factor <= 1.1d0 ) cell_factor = 1.1d0
@@ -75,10 +77,7 @@ SUBROUTINE allocate_nlpot
   !
   nqxq = INT( ( (SQRT(ecutrho) + qnorm) / dq + 4) * cell_factor )
   lmaxq = 2*lmaxkb+1
-  !
   IF (lmaxq > 0) ALLOCATE (qrad( nqxq, nbetam*(nbetam+1)/2, lmaxq, nsp))
-  ALLOCATE (becsum( nhm * (nhm + 1)/2, nat, nspin))
-  if (tqr) ALLOCATE (ebecsum( nhm * (nhm + 1)/2, nat, nspin))
   !
   ! Calculate dimensions for array tab (including a possible factor
   ! coming from cell contraction during variable cell relaxation/MD)

PW/src/init_us_1.f90

@@ -99,15 +99,17 @@ subroutine init_us_1
        rot_ylm(n,n1)=CMPLX(1.0_dp/sqrt2,0.d0,kind=DP)
        rot_ylm(n,n1+1)=CMPLX(0.d0, 1.0_dp/sqrt2,kind=DP)
      enddo
-     fcoef=(0.d0,0.d0)
-     dvan_so = (0.d0,0.d0)
-     qq_so=(0.d0,0.d0)
-     qq_at  = 0.d0
-     qq_nt=0.d0
-  else
+  endif
+  if ( nhm > 0 ) then
+     if (lspinorb) then
+        fcoef=(0.d0,0.d0)
+        dvan_so = (0.d0,0.d0)
+        qq_so=(0.d0,0.d0)
+     else
+        dvan = 0.d0
+     end if
      qq_nt=0.d0
      qq_at  = 0.d0
-     dvan = 0.d0
   endif
   !
   !   For each pseudopotential we initialize the indices nhtol, nhtolm,
@@ -140,7 +142,7 @@ subroutine init_us_1
      !
      ! ijtoh map augmentation channel indexes ih and jh to composite
      ! "triangular" index ijh
-     ijtoh(:,:,nt) = -1
+     if ( nhm > 0 ) ijtoh(:,:,nt) = -1
      ijv = 0
      do ih = 1,nh(nt)
          do jh = ih,nh(nt)
@@ -301,9 +303,11 @@ subroutine init_us_1
   endif
 #endif
   ! finally we set the atomic specific qq_at matrices
-  do na=1, nat
-     qq_at(:,:, na) = qq_nt(:,:,ityp(na))
-  end do
+  if ( nhm > 0 ) then
+     do na=1, nat
+        qq_at(:,:, na) = qq_nt(:,:,ityp(na))
+     end do
+  end if
   !
   !     fill the interpolation table tab
   !

PW/src/newd.f90

@@ -238,11 +238,11 @@ SUBROUTINE newd( )
            !
         ELSE
            !
-           DO is = 1, nspin
-              !
-              deeq(1:nht,1:nht,na,is) = dvan(1:nht,1:nht,nt)
-              !
-           ENDDO
+           if ( nht > 0 ) THEN
+              DO is = 1, nspin
+                 deeq(1:nht,1:nht,na,is) = dvan(1:nht,1:nht,nt)
+              ENDDO
+           end if
            !
         ENDIF
         !

PW/src/sum_band.f90

@@ -60,8 +60,10 @@ SUBROUTINE sum_band()
   !
   CALL start_clock( 'sum_band' )
   !
-  becsum(:,:,:) = 0.D0
-  if (tqr) ebecsum(:,:,:) = 0.D0
+  if ( nhm > 0 ) then
+     becsum(:,:,:) = 0.D0
+     if (tqr) ebecsum(:,:,:) = 0.D0
+  end if
   rho%of_r(:,:)      = 0.D0
   rho%of_g(:,:)      = (0.D0, 0.D0)
   if ( dft_is_meta() .OR. lxdm ) then