diff --git a/PW/cgramg1.f90 b/PW/cgramg1.f90
index c4416631a..092839358 100644
--- a/PW/cgramg1.f90
+++ b/PW/cgramg1.f90
@@ -89,21 +89,25 @@ SUBROUTINE cgramg1( lda, nvecx, n, start, finish, psi, spsi, hpsi )
        !
        DO vec = start, finish
           !
-	  CALL DGEMV( 'T', 2*n, vec-1, 2.D0, &
-	              psi, 2*n, spsi(1,vec), 1, 0.D0, ps, 1 )
-          !
-          IF ( gstart == 2 ) &
-	     ps(1:vec-1) = ps(1:vec-1) - psi(1,1:vec-1) * spsi(1,vec)
-          !
-          CALL reduce( ( vec - 1 ), ps )
-	  !
-          DO vecp = 1, ( vec - 1 )
+	  IF ( vec > 1 ) THEN
+	     !
+	     CALL DGEMV( 'T', 2*n, vec-1, 2.D0, &
+	                 psi, 2*n, spsi(1,vec), 1, 0.D0, ps, 1 )
              !
-             psi(:,vec)  = psi(:,vec)  - ps(vecp) * psi(:,vecp)
-             hpsi(:,vec) = hpsi(:,vec) - ps(vecp) * hpsi(:,vecp)
-             spsi(:,vec) = spsi(:,vec) - ps(vecp) * spsi(:,vecp)
+             IF ( gstart == 2 ) &
+	        ps(1:vec-1) = ps(1:vec-1) - psi(1,1:vec-1) * spsi(1,vec)
              !
-          END DO
+             CALL reduce( ( vec - 1 ), ps )
+	     !
+             DO vecp = 1, ( vec - 1 )
+                !
+                psi(:,vec)  = psi(:,vec)  - ps(vecp) * psi(:,vecp)
+                hpsi(:,vec) = hpsi(:,vec) - ps(vecp) * hpsi(:,vecp)
+                spsi(:,vec) = spsi(:,vec) - ps(vecp) * spsi(:,vecp)
+                !
+             END DO
+	     !
+	  END IF
           !
           psi_norm = 2.D0 * DDOT( 2 * n, psi(1,vec), 1, spsi(1,vec), 1 )
           !
@@ -160,17 +164,22 @@ SUBROUTINE cgramg1( lda, nvecx, n, start, finish, psi, spsi, hpsi )
        !
        DO vec = start, finish
           !
-	  CALL ZGEMV( 'C', n, vec-1, ONE, psi, n, spsi(1,vec), 1, ZERO, ps, 1 )
-	  !
-          CALL reduce( 2*( vec - 1 ), ps )
-          !
-          DO vecp = 1, ( vec - 1 )
+	  IF ( vec > 1 ) THEN
+	     !
+	     CALL ZGEMV( 'C', n, vec-1, ONE, &
+	                 psi, n, spsi(1,vec), 1, ZERO, ps, 1 )
+	     !
+             CALL reduce( 2*( vec - 1 ), ps )
              !
-             psi(:,vec)  = psi(:,vec)  - ps(vecp) * psi(:,vecp)
-             hpsi(:,vec) = hpsi(:,vec) - ps(vecp) * hpsi(:,vecp)
-             spsi(:,vec) = spsi(:,vec) - ps(vecp) * spsi(:,vecp)
-             !
-          END DO
+             DO vecp = 1, ( vec - 1 )
+                !
+                psi(:,vec)  = psi(:,vec)  - ps(vecp) * psi(:,vecp)
+                hpsi(:,vec) = hpsi(:,vec) - ps(vecp) * hpsi(:,vecp)
+                spsi(:,vec) = spsi(:,vec) - ps(vecp) * spsi(:,vecp)
+                !
+             END DO
+	     !
+	  END IF
           !
           psi_norm = DDOT( 2*n, psi(1,vec), 1, spsi(1,vec), 1 )
           !