change to pcg.f90 in ParO to force the residual and search vectors to be...

change to pcg.f90 in ParO to force the residual and search vectors to be orthogonal to the occupied manyfold.
Not totally sure of the S|psi0Xpsi0| or |psi0Xpsi0|S order

KS_Solvers/ParO/pcg.f90

@@ -50,6 +50,7 @@ SUBROUTINE pcg( hs_1psi, g_1psi, psi0, spsi0, npw, npwx, nbnd, psi, ethr, iter,
   !  ... solve the linear system
   !
   !      [ H - e S + lambda Pv ]|\tilde\psi> = [e S - H ] |psi>
+  !      Pc [ H - e S ]|\tilde\psi> = Pc [ e S - H ] |psi>
   !
   ! the solution is sought until the residual norm is a fixed fraction of the RHS norm
   ! in this way the more accurate is the original problem the more accuratly the correction is computed
@@ -154,6 +155,13 @@ SUBROUTINE pcg( hs_1psi, g_1psi, psi0, spsi0, npw, npwx, nbnd, psi, ethr, iter,
 
      psi(:) = psi(:) + alpha * p(:)
      r(:) = r(:) - alpha * w(:)
+!- project on conduction bands
+  CALL start_clock( 'pcg:ortho' )
+  CALL DGEMV( 'T', npw2, nbnd, 2.0_DP, psi0, npwx2,       r, 1, 0.0_DP, spsi0vec, 1 )
+  IF ( gstart == 2 ) spsi0vec(:) = spsi0vec(:) - CONJG(psi0(1,:))*r(1)
+  CALL mp_sum( spsi0vec, intra_bgrp_comm )
+  CALL DGEMV( 'N', npw2, nbnd, -1.D0, spsi0, npwx2, spsi0vec, 1, 1.0_DP,       r, 1 )
+  CALL stop_clock( 'pcg:ortho' )
      z(:) = r(:) ; call g_1psi(npwx,npw,z,e)
 !- project on conduction bands
   CALL start_clock( 'pcg:ortho' )