- FPMD/CP ortho subroutines merged as much as possible.

- Many ortho auxiliary functions (tauset, rhoset, sigset, calphi, updatc)
  are now in common between FPMD/CP, and moved to module ortho_base.f90
- In FPMD, three index vectors, related to real space like charge and potential
  have been substituted with single index vector like in CP, for compatibility
  and efficiency.
- Bug fix in pwtools/matdyn.f90 a logical variable was used in place of a
  character variable


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2694 c92efa57-630b-4861-b058-cf58834340f0

CPV/cg_sub.f90

@@ -25,7 +25,7 @@
       use energies, only: eht, epseu, exc, etot, eself, enl, ekin,          &
      &                    atot, entropy, egrand
       use electrons_base, only: f, nspin, nel, iupdwn, nupdwn, nudx, nelt, &
-                                nx => nbspx, n => nbsp, ispin => fspin
+                                nx => nbspx, n => nbsp, ispin
 
       use ensemble_dft, only: tens, tgrand, ninner, ismear, etemp, ef,       &
      &                tdynz, tdynf, zmass, fmass, fricz, fricf, z0, c0diag,  &
@@ -74,6 +74,7 @@
             berry_energy, ctabin, gqq, gqqm, df, pberryel
       use mp, only: mp_sum,mp_bcast
       use cp_electronic_mass,       ONLY : emass_cutoff
+      use orthogonalize_base,       ONLY : calphi
 
 !
       implicit none
@@ -143,7 +144,8 @@
 
       !calculates phi for pcdaga
 
-      call calphiid(c0,bec,betae,phi)
+      ! call calphiid(c0,bec,betae,phi)
+      CALL calphi( c0, SIZE(c0,1), bec, nhsa, betae, phi, n )
 
       !calculates the factors for S and K inversion in US case
       if(nvb.gt.0) then
@@ -694,7 +696,8 @@
  
         call calbec (1,nsp,eigr,c0,bec)
         !calculates phi for pc_daga
-        call calphiid(c0,bec,betae,phi)
+        !call calphiid(c0,bec,betae,phi)
+        CALL calphi( c0, SIZE(c0,1), bec, nhsa, betae, phi, n )
   
         !=======================================================================
         !

CPV/cglib.f90

@@ -6,93 +6,6 @@
 ! or http://www.gnu.org/copyleft/gpl.txt .
 !
 #include "f_defs.h"
-!-------------------------------------------------------------------------
-      subroutine calphiid(c0,bec,betae,phi)
-!-----------------------------------------------------------------------
-!     input: c0 (orthonormal with s(r(t)), bec=<c0|beta>, betae=|beta>
-!     computes the matrix phi (with the old positions)
-!       where  |phi> = s'|c0> = |c0> + sum q_ij |i><j|c0>
-!     where s'=s(r(t))  
-!     
-      use ions_base, only: na, nsp
-      use io_global, only: stdout
-      use cvan
-      use uspp_param, only: nh
-      use uspp, only :nhsa=>nkb, nhsavb=>nkbus, qq
-      use electrons_base, only: n => nbsp
-      use gvecw, only: ngw
-      use constants, only: pi, fpi
-      use control_flags, only: iprint, iprsta
-      use mp, only: mp_sum
-!
-      implicit none
-      complex(8) c0(ngw,n), phi(ngw,n), betae(ngw,nhsa)
-      real(8)    bec(nhsa,n) ,emtot
-! local variables
-      integer is, iv, jv, ia, inl, jnl, i, j
-      real(8) qtemp(nhsavb,n) ! automatic array
-!
-      phi(1:ngw,1:n) = 0.0d0
-!
-      if (nvb.gt.0) then
-         qtemp = 0.0d0
-         do is=1,nvb
-            do iv=1,nh(is)
-               do jv=1,nh(is)
-                  if(abs(qq(iv,jv,is)).gt.1.e-5) then
-                     do ia=1,na(is)
-                        inl=ish(is)+(iv-1)*na(is)+ia
-                        jnl=ish(is)+(jv-1)*na(is)+ia
-                        do i=1,n
-                           qtemp(inl,i) = qtemp(inl,i) +                &
-     &                                    qq(iv,jv,is)*bec(jnl,i)
-                        end do
-                     end do
-                  endif
-               end do
-            end do
-         end do
-         !
-         call MXMA                                                     &
-     &       (betae,1,2*ngw,qtemp,1,nhsavb,phi,1,2*ngw,2*ngw,nhsavb,n)
-      end if
-!
-      do j=1,n
-         do i=1,ngw
-            phi(i,j)=(phi(i,j)+c0(i,j))
-         end do
-      end do
-!     =================================================================
-      if(iprsta.gt.2) then
-         emtot=0.
-         do j=1,n
-            do i=1,ngw
-               emtot=emtot                                              &
-     &        +2.*DBLE(phi(i,j)*CONJG(c0(i,j)))
-            end do
-         end do
-         emtot=emtot/n
-         call mp_sum(emtot)
-         WRITE( stdout,*) 'in calphi sqrt(emtot)=',sqrt(emtot)
-         WRITE( stdout,*)
-         do is=1,nsp
-            if(nsp.gt.1) then
-               WRITE( stdout,'(33x,a,i4)') ' calphi: bec (is)',is
-               WRITE( stdout,'(8f9.4)')                                       &
-     &            ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)
-            else
-               do ia=1,na(is)
-                  WRITE( stdout,'(33x,a,i4)') ' calphi: bec (ia)',ia
-                  WRITE( stdout,'(8f9.4)')                                    &
-     &               ((bec(ish(is)+(iv-1)*na(is)+ia,i),iv=1,nh(is)),i=1,n)
-               end do
-            end if
-         end do
-      endif
-!
-      return
-    end subroutine calphiid
-
 
 
 !-----------------------------------------------------------------------
@@ -176,44 +89,26 @@
 !
       use parallel_toolkit, only: dspev_drv
       implicit none
-      integer nx,n,naux,ndim,iopt,iflag,k,i,j,info
+      integer nx,n,ndim,iflag,k,i,j
       real(8)   dval(n)
       real(8) amat(nx,n), dvec(nx,n)
       real(8), allocatable::  ap(:)
-      real(8), allocatable::  aux(:)
-
 
       ndim=(n*(n+1))/2
-      naux=3*n
       allocate(ap(ndim))
-      allocate(aux(naux))
       ap(:)=0.d0
-      aux(:)=0.d0
-
-      if (iflag.eq.1) then
-       iopt=1
-      else if (iflag.eq.0) then
-       iopt=0
-      else
-       write(*,*) 'ERROR: diag, iflag not allowed'
-       stop
-      end if
 
       k=0
       do j=1,n
        do i=1,j
         k=k+1
-!       ap(i + (j-1)*j/2)=amat(i,j)
         ap(k)=amat(i,j)
        end do
       end do
 
       CALL dspev_drv( 'V', 'U', n, ap, dval, dvec, nx )
 
-      if(info.ne.0) write(6,*) 'Problems with ddiag'
-
       deallocate(ap)
-      deallocate(aux)
 
       return
     end subroutine ddiag
@@ -293,7 +188,7 @@ subroutine pc2(a,beca,b,becb)
       use control_flags, only: iprint, iprsta
       use reciprocal_vectors, only: ng0 => gstart
       use mp, only: mp_sum
-      use electrons_base, only: n => nbsp, fspin
+      use electrons_base, only: n => nbsp, ispin
       use uspp_param, only: nh
       use uspp, only :nhsa=>nkb
       use uspp, only :qq
@@ -312,7 +207,7 @@ subroutine pc2(a,beca,b,becb)
             if (ng0.eq.2) then
                b(1,i)=0.5d0*(b(1,i)+CONJG(b(1,i)))
             endif
-            if(fspin(j) == fspin(i)) then
+            if(ispin(j) == ispin(i)) then
                do  ig=1,ngw           !loop on g vectors
                   sca=sca+2.d0*DBLE(CONJG(a(ig,j))*b(ig,i)) !2. for real wavefunctions
                enddo
@@ -369,7 +264,7 @@ subroutine pc2(a,beca,b,becb)
       use control_flags, only: iprint, iprsta
       use reciprocal_vectors, only: ng0 => gstart
       use mp, only: mp_sum
-      use electrons_base, only: n => nbsp, fspin
+      use electrons_base, only: n => nbsp, ispin
       use uspp_param, only: nh
       use uspp, only :nhsa=>nkb
 
@@ -387,7 +282,7 @@ subroutine pc2(a,beca,b,becb)
             if (ng0.eq.2) then
                b(1,i)=0.5d0*(b(1,i)+CONJG(b(1,i)))
             endif
-            if(fspin(i) == fspin(j)) then
+            if(ispin(i) == ispin(j)) then
                do  ig=1,ngw           !loop on g vectors
                   sca=sca+2.*DBLE(CONJG(a(ig,j))*b(ig,i)) !2. for real weavefunctions
                enddo
@@ -422,7 +317,7 @@ subroutine pc2(a,beca,b,becb)
       use control_flags, only: iprint, iprsta
       use reciprocal_vectors, only: ng0 => gstart
       use mp, only: mp_sum, mp_bcast
-      use electrons_base, only: n => nbsp, fspin
+      use electrons_base, only: n => nbsp, ispin
       use uspp_param, only: nh
       use uspp, only :nhsa=>nkb,qq, nhsavb=>nkbus
       use io_global, ONLY: ionode, ionode_id
@@ -535,7 +430,7 @@ subroutine pc2(a,beca,b,becb)
       use control_flags, only: iprint, iprsta
       use reciprocal_vectors, only: ng0 => gstart
       use mp, only: mp_sum, mp_bcast
-      use electrons_base, only: n => nbsp, fspin
+      use electrons_base, only: n => nbsp, ispin
       use uspp_param, only: nh
       use uspp, only :nhsa=>nkb,qq,nhsavb=>nkbus
       use io_global, ONLY: ionode, ionode_id

CPV/chargedensity.f90

@@ -32,7 +32,7 @@
 !  end of module-scope declarations
 !  ----------------------------------------------
 
-        PUBLIC :: checkrho, rhoofr, gradrho
+        PUBLIC :: rhoofr, gradrho
 
 !=----------------------------------------------------------------------=!
       CONTAINS
@@ -44,73 +44,6 @@
           RETURN
         END SUBROUTINE charge_density_closeup
 !
-!=----------------------------------------------------------------------=!
-    SUBROUTINE checkrho(rhoe, desc, rsum, omega)
-!=----------------------------------------------------------------------=!
-
-!     This Subroutine checks the value of the charge density integral
-!     that should be equal to the total charge 
-
-      USE constants, ONLY: rhothr
-      USE mp_global, ONLY: group, root, mpime
-      USE io_global, ONLY: ionode, stdout
-      USE mp, ONLY: mp_sum
-      USE charge_types, ONLY: charge_descriptor
-
-
-        IMPLICIT NONE
-
-        REAL(DP), INTENT(IN) :: omega 
-        REAL(DP) :: rsum(:)
-        REAL(DP), INTENT(IN) :: rhoe(:,:,:,:)
-        TYPE (charge_descriptor), INTENT(IN) :: desc
-        REAL(DP) :: rsum1
-        INTEGER :: i, j, k, ispin, nspin, nr1, nr2, nr3, ierr
-        INTEGER :: nxl, nyl, nzl
-
-        nr1 = desc%nx
-        nr2 = desc%ny
-        nr3 = desc%nz
-        nxl = desc%nxl
-        nyl = desc%nyl
-        nzl = desc%nzl
-        nspin  = desc%nspin
-
-! ...   recompute the integral of the charge density (for checking purpose)
-
-        DO ispin = 1, nspin
-
-          rsum1 = SUM( rhoe( 1:nxl, 1:nyl, 1:nzl, ispin ) )
-          rsum1 = rsum1 * omega / DBLE( nr1 * nr2 * nr3  )
-
-! ...     sum over all processors
-
-          CALL mp_sum( rsum1, group )
-          CALL mp_sum( rsum(ispin), group )
-
-! ...     write result (only processor 0)
-
-          IF( ionode ) THEN
-
-            WRITE( stdout,1) rsum(ispin), rsum1
-
-! ...       issue a warning if the result has changed
-
-            IF( ABS( rsum(ispin) - rsum1 ) > rhothr ) WRITE( stdout,100)
-
-          END IF
-
-        END DO
-
-    1 FORMAT(//,3X,'Total integrated electronic density',/ &
-     &         ,3X,'in G-space =',F11.6,4X,'in R-space =',F11.6)
-  100 FORMAT('** WARNING: CHARGE DENSITY **')
-
-        RETURN
-!=----------------------------------------------------------------------=!
-      END SUBROUTINE checkrho
-!=----------------------------------------------------------------------=!
-
 
       REAL(DP) FUNCTION dft_total_charge( ispin, c, cdesc, fi )
 
@@ -166,7 +99,7 @@
 !=----------------------------------------------------------------------=!
 !  BEGIN manual
 
-   SUBROUTINE rhoofr (nfi, c0, cdesc, fi, rhoe, desc, box)
+   SUBROUTINE rhoofr (nfi, c0, cdesc, fi, rhoe, box)
 
 !  this routine computes:
 !  rhoe = normalized electron density in real space
@@ -198,11 +131,12 @@
     USE turbo,           ONLY: tturbo, nturbo, turbo_states, allocate_turbo
     USE cell_module,     ONLY: boxdimensions
     USE wave_types,      ONLY: wave_descriptor
-    USE charge_types, ONLY: charge_descriptor
     USE io_global,       ONLY: stdout, ionode
     USE control_flags,   ONLY: force_pairing, iprint
     USE parameters,      ONLY: nspinx
     USE brillouin,       ONLY: kpoints, kp
+    USE grid_dimensions, ONLY: nr1, nr2, nr3, nr1x, nr2x, nnrx
+
 
 
     IMPLICIT NONE
@@ -213,32 +147,23 @@
     COMPLEX(DP)                     :: c0(:,:,:,:)
     TYPE (boxdimensions), INTENT(IN) :: box
     REAL(DP),          INTENT(IN) :: fi(:,:,:)
-    REAL(DP),            INTENT(OUT) :: rhoe(:,:,:,:)
+    REAL(DP),            INTENT(OUT) :: rhoe(:,:)
-    TYPE (charge_descriptor),    INTENT(IN) :: desc
     TYPE (wave_descriptor), INTENT(IN) :: cdesc
 
 ! ... declare other variables
 
-    INTEGER :: i, is1, is2, j, k, ib, ik, nb, nxl, nyl, nzl, ispin
+    INTEGER :: i, is1, is2, j, k, ib, ik, nb, ispin
-    INTEGER :: nr1x, nr2x, nr3x, nspin, nbnd, nnr
+    INTEGER :: nspin, nbnd, nnr
     REAL(DP)  :: r2, r1, coef3, coef4, omega, rsumg( nspinx ), rsumgs
     REAL(DP)  :: fact, rsumr( nspinx )
-    REAL(DP), ALLOCATABLE :: rho(:,:,:)
+    COMPLEX(DP), ALLOCATABLE :: psi2(:)
-    COMPLEX(DP), ALLOCATABLE :: psi2(:,:,:)
     INTEGER :: ierr, ispin_wfc
     LOGICAL :: ttprint
 
 ! ... end of declarations
 !  ----------------------------------------------
 
-    nxl =  dfftp%nr1
+    nnr   =  dfftp%nr1x * dfftp%nr2x * dfftp%npl
-    nyl =  dfftp%nr2
-    nzl =  dfftp%npl
-    nnr =  dfftp%nr1 * dfftp%nr2 * dfftp%nr3
-
-    nr1x = dfftp%nr1x
-    nr2x = dfftp%nr2x
-    nr3x = dfftp%npl
 
     omega = box%deth
 
@@ -250,29 +175,20 @@
     ttprint = .FALSE.
     IF( nfi == 0 .or. mod( nfi, iprint ) == 0 ) ttprint = .TRUE.
 
-    ! ... Check consistensy of the charge density grid and fft grid
+    ALLOCATE( psi2( nnrx ), STAT=ierr )
-
-    IF( SIZE( rhoe, 1 ) < nxl ) &
-      CALL errore(' rhoofr ', ' wrong X dimension for rhoe ',1)
-    IF( SIZE( rhoe, 2 ) < nyl ) &
-      CALL errore(' rhoofr ', ' wrong Y dimension for rhoe ',1)
-    IF( SIZE( rhoe, 3 ) < nzl ) &
-      CALL errore(' rhoofr ', ' wrong Z dimension for rhoe ',1)
-
-    ALLOCATE( psi2( nr1x, nr2x, nr3x ), STAT=ierr )
     IF( ierr /= 0 ) CALL errore(' rhoofr ', ' allocating psi2 ', ABS(ierr) )
-    ALLOCATE( rho( nr1x, nr2x, nr3x ), STAT=ierr )
-    IF( ierr /= 0 ) CALL errore(' rhoofr ', ' allocating rho ', ABS(ierr) )
 
     IF( tturbo ) THEN
       !
       ! ... if tturbo=.TRUE. some data is stored in memory instead of being
       ! ... recalculated (see card 'TURBO')
       !
-      CALL allocate_turbo( dfftp%nr1x, dfftp%nr2x, dfftp%npl )
+      CALL allocate_turbo( nnrx )
 
     END IF
 
+    rhoe  = zero
+
     DO ispin = 1, nspin
 
       ! ... arrange for FFT of wave functions
@@ -285,7 +201,7 @@
         ! ...  Gamma-point calculation: wave functions are real and can be
         ! ...  Fourier-transformed two at a time as a complex vector
 
-        rho = zero
+        psi2 = zero
 
         nbnd = cdesc%nbl( ispin )
         nb   = ( nbnd - MOD( nbnd, 2 ) )
@@ -302,7 +218,7 @@
 
           IF( tturbo .AND. ( ib <= nturbo ) ) THEN
             ! ...  store real-space wave functions to be used in force 
-            turbo_states( :, :, :, ib ) = psi2( :, :, : )
+            turbo_states( :, ib ) = psi2( : )
           END IF
 
           ! ...  occupation numbers divided by cell volume
@@ -313,22 +229,18 @@
 
           ! ...  compute charge density from wave functions
 
-          DO k = 1, nzl
+          DO i = 1, nnr
-            DO j = 1, nyl
-              DO i = 1, nxl
 
                 ! ...  extract wave functions from psi2
 
-                r1 =  DBLE( psi2(i,j,k) ) 
+                r1 =  DBLE( psi2(i) ) 
-                r2 = AIMAG( psi2(i,j,k) ) 
+                r2 = AIMAG( psi2(i) ) 
 
                 ! ...  add squared moduli to charge density
 
-                rho(i,j,k) = rho(i,j,k) + coef3 * r1 * r1 + coef4 * r2 * r2
+                rhoe(i,ispin) = rhoe(i,ispin) + coef3 * r1 * r1 + coef4 * r2 * r2
 
           END DO
-            END DO
-          END DO
 
         END DO
 
@@ -346,21 +258,17 @@
 
           ! ...  compute charge density from wave functions
 
-          DO k = 1, nzl
+          DO i = 1, nnr
-            DO j = 1, nyl
-              DO i = 1, nxl
 
              ! ...  extract wave functions from psi2
 
-                r1 = DBLE( psi2(i,j,k) )
+             r1 = DBLE( psi2(i) )
 
              ! ...  add squared moduli to charge density
 
-                rho(i,j,k) = rho(i,j,k) + coef3 * r1 * r1
+             rhoe(i,ispin) = rhoe(i,ispin) + coef3 * r1 * r1
 
           END DO
-            END DO
-          END DO
 
         END IF
 
@@ -368,7 +276,7 @@
 
         ! ...  calculation with generic k points: wave functions are complex
 
-        rho = zero
+        psi2 = zero
 
         DO ik = 1, cdesc%nkl
 
@@ -385,25 +293,19 @@
 
             ! ...  compute charge density
 
-            DO k = 1, nzl
+            DO i = 1, nnr
-              DO j = 1, nyl
-                DO i = 1, nxl
 
                 ! ...  add squared modulus to charge density
 
-                  rho(i,j,k) = rho(i,j,k) + coef3 * DBLE( psi2(i,j,k) * CONJG(psi2(i,j,k)) )
+                rhoe(i,ispin) = rhoe(i,ispin) + coef3 * DBLE( psi2(i) * CONJG(psi2(i)) )
 
             END DO
           END DO
         END DO
-          END DO
-        END DO
 
       END IF
 
-      IF( ttprint ) rsumr( ispin ) = SUM( rho ) * omega / nnr
+      IF( ttprint ) rsumr( ispin ) = SUM( rhoe( :, ispin ) ) * omega / ( nr1 * nr2 * nr3 )
-
-      rhoe( 1:nxl, 1:nyl, 1:nzl, ispin ) = rho( 1:nxl, 1:nyl, 1:nzl )
 
     END DO
 
@@ -443,8 +345,6 @@
 
     DEALLOCATE(psi2, STAT=ierr)
     IF( ierr /= 0 ) CALL errore(' rhoofr ', ' deallocating psi2 ', ABS(ierr) )
-    DEALLOCATE(rho, STAT=ierr)
-    IF( ierr /= 0 ) CALL errore(' rhoofr ', ' deallocating rho ', ABS(ierr) )
 
 
     RETURN
@@ -468,7 +368,7 @@
 
      COMPLEX(DP), INTENT(IN)  :: rhoeg(:)    ! charge density (Reciprocal Space)
      REAL(DP), INTENT(IN)  :: gx(:,:)        ! cartesian components of G-vectors
-     REAL(DP),  INTENT(OUT) :: grho(:,:,:,:) ! charge density gradient
+     REAL(DP),  INTENT(OUT) :: grho(:,:) ! charge density gradient
 
      INTEGER :: ig, ipol, ierr
      COMPLEX(DP), ALLOCATABLE :: tgrho(:)
@@ -484,7 +384,7 @@
          rg        = rhoeg(ig) * gx( ipol, ig )
          tgrho(ig) = tpiba * CMPLX( - AIMAG(rg), DBLE(rg) ) 
        END DO
-       CALL pinvfft( grho(:,:,:,ipol), tgrho )
+       CALL pinvfft( grho(:,ipol), tgrho )
      END DO
 
      DEALLOCATE(tgrho, STAT=ierr)

CPV/chargemix.f90

@@ -160,7 +160,7 @@
       IMPLICIT NONE
 
 ! ... declare subroutine arguments
-      REAL(DP), INTENT(INOUT) :: rhoe(:,:,:)
+      REAL(DP), INTENT(INOUT) :: rhoe(:)
       REAL(DP), INTENT(OUT) ::  drho
       INTEGER, INTENT(IN) :: nfi
 
@@ -170,7 +170,7 @@
       REAL(DP) :: g02, g12, ar, den, num, rsc
       REAL(DP) :: alpha(daamax)
       REAL(DP), ALLOCATABLE :: aa(:,:)
-      REAL(DP), ALLOCATABLE :: rho_old(:,:,:)
+      REAL(DP), ALLOCATABLE :: rho_old(:)
       INTEGER :: ns, sp, is, ism, i, ig
       LOGICAL, SAVE :: tfirst = .TRUE.
       INTEGER, SAVE :: dimaa, dimaaold, nrho_t, ierr
@@ -305,7 +305,7 @@
 
       END IF
 
-      ALLOCATE( rho_old( SIZE(rhoe, 1), SIZE(rhoe, 2), SIZE(rhoe, 3) ), STAT=ierr )
+      ALLOCATE( rho_old( SIZE(rhoe) ), STAT=ierr )
       IF( ierr /= 0 ) CALL errore(' newrho ', ' allocating rho_old ', ierr)
       rho_old = rhoe 
 

CPV/cp_restart.f90

@@ -50,7 +50,7 @@ MODULE cp_restart
       USE gvecw,                    ONLY : ngw, ngwt, ecutw, gcutw
       USE reciprocal_vectors,       ONLY : ig_l2g, mill_l
       USE electrons_base,           ONLY : nspin, nbnd, nbsp, nelt, nel, &
-                                           nupdwn, iupdwn, f, fspin, nudx
+                                           nupdwn, iupdwn, f, nudx
       USE cell_base,                ONLY : ibrav, alat, celldm, &
                                            symm_type, s_to_r
       USE ions_base,                ONLY : nsp, nat, na, atm, zv, &
@@ -1129,6 +1129,7 @@ MODULE cp_restart
          IF ( ( nspin_ /= nspin ) .OR. &
               ( nbnd_ /= nbnd ) .OR. ( NINT( nelec_ ) /= nelt ) ) THEN 
             !
+            attr = "electron, bands or spin do not match"
             ierr = 30
             !
             GOTO 100

CPV/cplib.f90

@@ -238,108 +238,6 @@
 !
 
 
-!-------------------------------------------------------------------------
-      SUBROUTINE calphi( c0, ngwx, ema0bg, bec, nkbx, betae, phi, n )
-!-----------------------------------------------------------------------
-!     input: c0 (orthonormal with s(r(t)), bec=<c0|beta>, betae=|beta>
-!     computes the matrix phi (with the old positions)
-!       where  |phi> = s'|c0> = |c0> + sum q_ij |i><j|c0>
-!     where s'=s(r(t))  
-!
-      USE kinds,          ONLY: DP
-      USE ions_base,      ONLY: na, nsp
-      USE io_global,      ONLY: stdout
-      USE cvan,           ONLY: ish, nvb
-      USE uspp_param,     ONLY: nh
-      USE uspp,           ONLY: nhsavb=>nkbus, qq
-      USE gvecw,          ONLY: ngw
-      USE constants,      ONLY: pi, fpi
-      USE control_flags,  ONLY: iprint, iprsta
-      USE mp,             ONLY: mp_sum
-!
-      IMPLICIT NONE
-      
-      INTEGER, INTENT(IN) :: ngwx, nkbx, n
-      COMPLEX(DP)         :: c0( ngwx, n ), phi( ngwx, n ), betae( ngwx, nkbx )
-      REAL(DP)            :: ema0bg( ngwx ), bec( nkbx, n ), emtot
-
-      ! local variables
-      !
-      INTEGER  :: is, iv, jv, ia, inl, jnl, i, j
-      REAL(DP), ALLOCATABLE :: qtemp( : , : )
-!
-      CALL start_clock( 'calphi' )
-
-      ALLOCATE( qtemp( nhsavb, n ) )
-
-      phi(:,:) = (0.d0, 0.d0)
-!
-      IF ( nvb > 0 ) THEN
-         qtemp (:,:) = 0.d0
-         DO is=1,nvb
-            DO iv=1,nh(is)
-               DO jv=1,nh(is)
-                  IF(ABS(qq(iv,jv,is)) > 1.e-5) THEN
-                     DO ia=1,na(is)
-                        inl=ish(is)+(iv-1)*na(is)+ia
-                        jnl=ish(is)+(jv-1)*na(is)+ia
-                        DO i=1,n
-                           qtemp(inl,i) = qtemp(inl,i) +                &
-     &                                    qq(iv,jv,is)*bec(jnl,i)
-                        END DO
-                     END DO
-                  ENDIF
-               END DO
-            END DO
-         END DO
-!
-         CALL MXMA                                                     &
-     &       ( betae, 1, 2*ngwx, qtemp, 1, nhsavb, phi, 1, 2*ngwx, 2*ngw, nhsavb, n )
-      END IF
-!
-      DO j=1,n
-         DO i=1,ngw
-            phi(i,j)=(phi(i,j)+c0(i,j))*ema0bg(i)
-         END DO
-      END DO
-!     =================================================================
-      IF(iprsta > 2) THEN
-         emtot=0.0d0
-         DO j=1,n
-            DO i=1,ngw
-               emtot=emtot                                              &
-     &        +2.0d0*DBLE(phi(i,j)*CONJG(c0(i,j)))*ema0bg(i)**(-2.0d0)
-            END DO
-         END DO
-         emtot=emtot/n
-
-         CALL mp_sum( emtot )
-
-         WRITE( stdout,*) 'in calphi sqrt(emtot)=',SQRT(emtot)
-         WRITE( stdout,*)
-         DO is=1,nsp
-            IF(nsp > 1) THEN
-               WRITE( stdout,'(33x,a,i4)') ' calphi: bec (is)',is
-               WRITE( stdout,'(8f9.4)')                                       &
-     &            ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)
-            ELSE
-               DO ia=1,na(is)
-                  WRITE( stdout,'(33x,a,i4)') ' calphi: bec (ia)',ia
-                  WRITE( stdout,'(8f9.4)')                                    &
-     &               ((bec(ish(is)+(iv-1)*na(is)+ia,i),iv=1,nh(is)),i=1,n)
-               END DO
-            END IF
-         END DO
-      ENDIF
-
-      DEALLOCATE( qtemp )
-
-      CALL stop_clock( 'calphi' )
-!
-      RETURN
-      END SUBROUTINE calphi
-
-
 
 !-----------------------------------------------------------------------
       REAL(8) FUNCTION cscnorm( bec, nkbx, cp, ngwx, i, n )
@@ -672,7 +570,7 @@
       USE uspp_param, ONLY: nhm, nh
       USE smooth_grid_dimensions, ONLY: nr1s, nr2s, nr3s, &
             nr1sx, nr2sx, nr3sx, nnrsx
-      USE electrons_base, ONLY: n => nbsp, ispin => fspin, f, nspin
+      USE electrons_base, ONLY: n => nbsp, ispin, f, nspin
       USE constants, ONLY: pi, fpi
       USE ions_base, ONLY: nsp, na, nat
       USE gvecw, ONLY: ggp
@@ -1151,26 +1049,34 @@
 
 !
 !-----------------------------------------------------------------------
-      REAL(8) FUNCTION enkin(c)
+   FUNCTION enkin( c, ngwx, f, n )
 !-----------------------------------------------------------------------
       !
       ! calculation of kinetic energy term
       !
+      USE kinds,              ONLY: DP
       USE constants,          ONLY: pi, fpi
-      USE electrons_base, ONLY: nx => nbspx, n => nbsp, f
       USE gvecw,              ONLY: ngw
       USE reciprocal_vectors, ONLY: gstart
       USE gvecw,              ONLY: ggp
       USE mp,                 ONLY: mp_sum
+      USE mp_global,          ONLY: group
       USE cell_base,          ONLY: tpiba2
 
       IMPLICIT NONE
+
+      REAL(DP)                :: enkin
+
       ! input
-      COMPLEX(8) c(ngw,nx)
+
-! local
+      INTEGER,     INTENT(IN) :: ngwx, n
-      INTEGER ig, i
+      COMPLEX(DP), INTENT(IN) :: c( ngwx, n )
-      REAL(8) sk(n)  ! automatic array
+      REAL(DP),    INTENT(IN) :: f( n )
       !
+      ! local
+
+      INTEGER  :: ig, i
+      REAL(DP) :: sk(n)  ! automatic array
       !
       DO i=1,n
          sk(i)=0.0
@@ -1179,12 +1085,16 @@
          END DO
       END DO
 
-      CALL mp_sum( sk(1:n) )
+      CALL mp_sum( sk(1:n), group )
 
       enkin=0.0
       DO i=1,n
          enkin=enkin+f(i)*sk(i)
       END DO
+
+      ! ... reciprocal-space vectors are in units of alat/(2 pi) so a
+      ! ... multiplicative factor (2 pi/alat)**2 is required
+
       enkin = enkin * tpiba2
 !
       RETURN
@@ -1464,7 +1374,7 @@
       USE cvan,           ONLY :nvb, ish
       USE uspp,           ONLY : nkb, nhsavb=>nkbus, qq
       USE uspp_param,     ONLY:  nh
-      USE electrons_base, ONLY: ispin => fspin
+      USE electrons_base, ONLY: ispin
       USE gvecw,          ONLY: ngw
       USE mp,             ONLY: mp_sum
       USE kinds,          ONLY: DP
@@ -2106,185 +2016,6 @@
       RETURN
       END SUBROUTINE nlfl
 
-
-!-----------------------------------------------------------------------
-      SUBROUTINE ortho                                                  &
-     &      (eigr,cp,phi,x0,diff,iter,ccc,eps,max,delt,bephi,becp)
-!-----------------------------------------------------------------------
-!     input = cp (non-orthonormal), beta
-!     input = phi |phi>=s'|c0>
-!     output= cp (orthonormal with s( r(t+dt) ) )
-!     output= bephi, becp
-!     the method used is similar to the version in les houches 1988
-!     'simple molecular systems at..'  p. 462-463  (18-22)
-!      xcx + b x + b^t x^t + a = 1
-!     where c = <s'c0|s|s'c0>   b = <s'c0|s cp>   a = <cp|s|cp>
-!     where s=s(r(t+dt)) and s'=s(r(t))  
-!     for vanderbilt pseudo pot - kl & ap
-!
-      USE ions_base, ONLY: na, nat
-      USE cvan, ONLY: ish, nvb
-      USE uspp, ONLY : nkb, qq
-      USE uspp_param, ONLY: nh
-      USE electrons_base, ONLY: n => nbsp, nbspx, nudx, nspin, nupdwn, iupdwn, f
-      USE gvecw, ONLY: ngw
-      USE control_flags, ONLY: iprint, iprsta
-      USE io_global, ONLY: stdout
-      USE orthogonalize_base, ONLY: ortho_iterate, diagonalize_rho, sigset, rhoset, &
-                                    tauset
-!
-      IMPLICIT NONE
-!
-      COMPLEX(8)   cp(ngw,n), phi(ngw,n), eigr(ngw,nat)
-      REAL(8) x0( nbspx, nbspx ), diff, ccc, eps, delt
-      INTEGER iter, max
-      REAL(8) bephi(nkb,n), becp(nkb,n)
-!
-      REAL(8), ALLOCATABLE :: diag(:), work1(:), work2(:), xloc(:,:), &
-                                   rhos(:,:), rhor(:,:), u(:,:), &
-                                   sig(:,:), rho(:,:), tau(:,:)
-
-      INTEGER :: ngwx, nkbx
-
-      INTEGER istart, nss, ifail, i, j, iss, iv, jv, ia, is, inl, jnl
-      REAL(8), ALLOCATABLE:: qbephi(:,:), qbecp(:,:)
-
-      ALLOCATE( diag( nudx ), work1( nudx ), work2( nudx ), xloc( nudx, nudx ),  &
-                rhos( nudx, nudx ), rhor( nudx, nudx ), u( nudx, nudx ),         &
-                sig( nudx, nudx ), rho( nudx, nudx ), tau( nudx, nudx ) )
-
-      ngwx = ngw
-      nkbx = nkb
-      !
-      !     calculation of becp and bephi
-      !
-      CALL start_clock( 'ortho' )
-
-      CALL nlsm1( n, 1, nvb, eigr,  cp,  becp )
-      CALL nlsm1( n, 1, nvb, eigr, phi, bephi )
-      !
-      !     calculation of qbephi and qbecp
-      !
-      ALLOCATE( qbephi( nkbx, n ) )
-      ALLOCATE( qbecp ( nkbx, n ) )
-      !
-      qbephi = 0.d0
-      qbecp  = 0.d0
-!
-      DO is=1,nvb
-         DO iv=1,nh(is)
-            DO jv=1,nh(is)
-               IF(ABS(qq(iv,jv,is)).GT.1.e-5) THEN
-                  DO ia=1,na(is)
-                     inl=ish(is)+(iv-1)*na(is)+ia
-                     jnl=ish(is)+(jv-1)*na(is)+ia
-                     DO i=1,n
-                        qbephi(inl,i)= qbephi(inl,i)                    &
-     &                       +qq(iv,jv,is)*bephi(jnl,i)
-                        qbecp (inl,i)=qbecp (inl,i)                     &
-     &                       +qq(iv,jv,is)*becp (jnl,i)
-                     END DO
-                  END DO
-               ENDIF
-            END DO
-         END DO
-      END DO
-!
-      DO iss = 1, nspin
-
-         nss    = nupdwn(iss)
-         istart = iupdwn(iss)
-
-         !
-         !     rho = <s'c0|s|cp>
-         !     sig = 1-<cp|s|cp>
-         !     tau = <s'c0|s|s'c0>
-         !
-         CALL rhoset( cp, ngwx, phi, bephi, nkbx, qbecp, n, nss, istart, rho, nudx )
-         !
-         CALL sigset( cp, ngwx, becp, nkbx, qbecp, n, nss, istart, sig, nudx )
-         !
-         CALL tauset( phi, ngwx, bephi, nkbx, qbephi, n, nss, istart, tau, nudx )
-         !
-!
-         IF(iprsta.GT.4) THEN
-            WRITE( stdout,*)
-            WRITE( stdout,'(26x,a)') '    rho '
-            DO i=1,nss
-               WRITE( stdout,'(7f11.6)') (rho(i,j),j=1,nss)
-            END DO
-            WRITE( stdout,*)
-            WRITE( stdout,'(26x,a)') '    sig '
-            DO i=1,nss
-               WRITE( stdout,'(7f11.6)') (sig(i,j),j=1,nss)
-            END DO
-            WRITE( stdout,*)
-            WRITE( stdout,'(26x,a)') '    tau '
-            DO i=1,nss
-               WRITE( stdout,'(7f11.6)') (tau(i,j),j=1,nss)
-            END DO
-         ENDIF
-!
-!
-!----------------------------------------------------------------by ap--
-! 
-         DO j=1,nss
-            DO i=1,nss
-               xloc(i,j) = x0(istart-1+i,istart-1+j)*ccc
-               rhos(i,j)=0.5d0*( rho(i,j)+rho(j,i) )
-!
-! on some machines (IBM RS/6000 for instance) the following test allows
-! to distinguish between Numbers and Sodium Nitride (NaN, Not a Number).
-! If a matrix of Not-Numbers is passed to rs, the most likely outcome is
-! that the program goes on forever doing nothing and writing nothing.
-!
-               IF (rhos(i,j).NE.rhos(i,j))                                &
-     &             CALL errore('ortho','ortho went bananas',1)
-               rhor(i,j)=rho(i,j)-rhos(i,j)
-            END DO
-         END DO
-!     
-         ifail=0
-         CALL start_clock( 'rsg' )
-         CALL diagonalize_rho( nss, rhos, diag, u )
-         ! CALL rs(nudx,nss,rhos,diag,1,u,work1,work2,ifail) 
-         CALL stop_clock( 'rsg' )
-!
-!                calculation of lagranges multipliers
-!
-         CALL ortho_iterate( u, diag, xloc, sig, rhor, rhos, tau, nudx, nss, max, eps )
-
-         IF(iprsta.GT.4) THEN
-            WRITE( stdout,*)
-            WRITE( stdout,'(26x,a)') '    lambda '
-            DO i=1,nss
-               WRITE( stdout,'(7f11.6)') (xloc(i,j)/f(i+istart-1),j=1,nss)
-            END DO
-         ENDIF
-!     
-         IF(iprsta.GT.2) THEN
-            WRITE( stdout,*) ' diff= ',diff,' iter= ',iter
-         ENDIF
-!     
-!     lagrange multipliers
-!
-         DO i=1,nss
-            DO j=1,nss
-               x0(istart-1+i,istart-1+j)=xloc(i,j)/ccc
-               IF (xloc(i,j).NE.xloc(i,j))                                &
-     &             CALL errore('ortho','ortho went bananas',2)
-            END DO
-         END DO
-!
-      END DO
-!
-      DEALLOCATE(qbecp )
-      DEALLOCATE(qbephi)
-      DEALLOCATE( diag, work1, work2, xloc, rhos, rhor, u, sig, rho, tau )
-!
-      CALL stop_clock( 'ortho' )
-      RETURN
-      END SUBROUTINE ortho
 !
 !-----------------------------------------------------------------------
       SUBROUTINE pbc(rin,a1,a2,a3,ainv,rout)
@@ -2624,6 +2355,7 @@
 !
       RETURN
       END SUBROUTINE rdiag
+
 !-----------------------------------------------------------------------
    SUBROUTINE rhoofr (nfi,c,irb,eigrb,bec,rhovan,rhor,rhog,rhos,enl,ekin)
 !-----------------------------------------------------------------------
@@ -2638,101 +2370,108 @@
 !
 !     e_v = sum_i,ij rho_i,ij d^ion_is,ji
 !
-      USE kinds, ONLY: dp
+      USE kinds,              ONLY: DP
       USE control_flags,      ONLY: iprint, tbuff, iprsta, thdyn, tpre, trhor
-      USE ions_base, ONLY: nat, nas => nax, nsp
+      USE ions_base,          ONLY: nat
-      USE parameters, ONLY: natx, nsx
+      USE gvecp,              ONLY: ngm
-      USE gvecp, ONLY: ng => ngm
+      USE gvecs,              ONLY: ngs, nps, nms
-      USE gvecs
       USE gvecb,              ONLY: ngb
       USE gvecw,              ONLY: ngw
-      USE reciprocal_vectors, ONLY: gstart
       USE recvecs_indexes,    ONLY: np, nm
-      USE uspp, ONLY: nhsa => nkb
+      USE reciprocal_vectors, ONLY: gstart
+      USE uspp,               ONLY: nkb
       USE uspp_param,         ONLY: nh, nhm
       USE grid_dimensions,    ONLY: nr1, nr2, nr3, &
-            nr1x, nr2x, nr3x, nnr => nnrx
+                                    nr1x, nr2x, nr3x, nnrx
       USE cell_base,          ONLY: omega
       USE smooth_grid_dimensions, ONLY: nr1s, nr2s, nr3s, &
                                         nr1sx, nr2sx, nr3sx, nnrsx
-      USE electrons_base, ONLY: nx => nbspx, n => nbsp, f, ispin => fspin, nspin
+      USE electrons_base,     ONLY: nx => nbspx, n => nbsp, f, ispin, nspin
       USE constants,          ONLY: pi, fpi
       USE mp,                 ONLY: mp_sum
-      ! use local_pseudo
+      USE dener,              ONLY: denl, dekin
-!
-      USE cdvan
-      USE dener
       USE io_global,          ONLY: stdout
       USE funct,              ONLY: dft_is_meta
       USE cg_module,          ONLY: tcg
+      USE cp_main_variables,  ONLY: rhopr
 !
       IMPLICIT NONE
-      REAL(8) bec(nhsa,n), rhovan(nhm*(nhm+1)/2,nat,nspin)
+      REAL(DP) bec(nkb,n), rhovan( nhm * ( nhm + 1 ) / 2, nat, nspin )
-      REAL(8) rhor(nnr,nspin), rhos(nnrsx,nspin)
+      REAL(DP) rhor(nnrx,nspin), rhos(nnrsx,nspin)
-      REAL(8) enl, ekin
+      REAL(DP) enl, ekin
-      COMPLEX(8) eigrb(ngb,nat), c(ngw,nx), rhog(ng,nspin)
+      COMPLEX(DP) eigrb( ngb, nat ), c( ngw, nx ), rhog( ngm, nspin )
       INTEGER irb( 3, nat ), nfi
+      
+
 ! local variables
+
       INTEGER iss, isup, isdw, iss1, iss2, ios, i, ir, ig
-      REAL(8) rsumr(2), rsumg(2), sa1, sa2
+      REAL(DP) rsumr(2), rsumg(2), sa1, sa2
-      REAL(8) rnegsum, rmin, rmax, rsum
+      REAL(DP) rnegsum, rmin, rmax, rsum
-      REAL(8), EXTERNAL :: enkin, ennl
+      REAL(DP), EXTERNAL :: enkin, ennl
-      COMPLEX(8) ci,fp,fm
+      COMPLEX(DP) ci,fp,fm
-      COMPLEX(8), ALLOCATABLE :: psi(:), psis(:)
+      COMPLEX(DP), ALLOCATABLE :: psi(:), psis(:)
-!
+      LOGICAL, SAVE :: first = .TRUE.
+      
       !
+
       CALL start_clock( 'rhoofr' )
-      ALLOCATE( psi( nnr ) ) 
+
+      ALLOCATE( psi( nnrx ) ) 
       ALLOCATE( psis( nnrsx ) ) 
+
       ci=(0.0,1.0)
+
       DO iss=1,nspin
          rhor(:,iss) = 0.d0
          rhos(:,iss) = 0.d0
          rhog(:,iss) = (0.d0, 0.d0)
       END DO
       !
-!     ==================================================================
       !  calculation of kinetic energy ekin
-!     ==================================================================
+      !
-      ekin=enkin(c)
+      ekin = enkin( c, ngw, f, n )
+      !
       IF( tpre ) CALL denkin( c, dekin )
       !
-!     ==================================================================
       !     calculation of non-local energy
-!     ==================================================================
+      !
       enl = ennl( rhovan, bec )
+      !
       IF( tpre ) CALL dennl( bec, denl )
       !    
       !    warning! trhor and thdyn are not compatible yet!   
       !
       IF( trhor .AND. ( .NOT. thdyn ) ) THEN
-!     ==================================================================
+         !
+         !   non self-consistent calculation  
          !   charge density is read from unit 47
-!     ==================================================================
+         !
-#ifdef __PARA
+         IF( first ) THEN
             CALL read_rho( 47, nspin, rhor )
-#else
+            rhopr = rhor
-         READ(47) ((rhor(ir,iss),ir=1,nnr),iss=1,nspin)
+            first = .FALSE.
-#endif
+         ELSE
-         REWIND 47
+            rhor = rhopr
+         END IF
 !
          IF(nspin.EQ.1)THEN
             iss=1
-            DO ir=1,nnr
+            DO ir=1,nnrx
                psi(ir)=CMPLX(rhor(ir,iss),0.d0)
             END DO
             CALL fwfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x)
-            DO ig=1,ng
+            DO ig=1,ngm
                rhog(ig,iss)=psi(np(ig))
             END DO
          ELSE
             isup=1
             isdw=2
-            DO ir=1,nnr
+            DO ir=1,nnrx
                psi(ir)=CMPLX(rhor(ir,isup),rhor(ir,isdw))
             END DO
             CALL fwfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x)
-            DO ig=1,ng
+            DO ig=1,ngm
                fp=psi(np(ig))+psi(nm(ig))
                fm=psi(np(ig))-psi(nm(ig))
                rhog(ig,isup)=0.5*CMPLX( DBLE(fp),AIMAG(fm))
@@ -2836,7 +2575,7 @@
                psi(np(ig))=      rhog(ig,iss)
             END DO
             CALL invfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x)
-            DO ir=1,nnr
+            DO ir=1,nnrx
                rhor(ir,iss)=DBLE(psi(ir))
             END DO
          ELSE 
@@ -2851,7 +2590,7 @@
                psi(np(ig))=rhog(ig,isup)+ci*rhog(ig,isdw)
             END DO
             CALL invfft(psi,nr1,nr2,nr3,nr1x,nr2x,nr3x)
-            DO ir=1,nnr
+            DO ir=1,nnrx
                rhor(ir,isup)= DBLE(psi(ir))
                rhor(ir,isdw)=AIMAG(psi(ir))
             END DO
@@ -2898,7 +2637,7 @@
 !
 !
       IF(iprsta.GE.2) THEN
-         CALL checkrho(nnr,nspin,rhor,rmin,rmax,rsum,rnegsum)
+         CALL checkrho(nnrx,nspin,rhor,rmin,rmax,rsum,rnegsum)
          rnegsum=rnegsum*omega/DBLE(nr1*nr2*nr3)
          rsum=rsum*omega/DBLE(nr1*nr2*nr3)
          WRITE( stdout,'(a,4(1x,f12.6))')                                     &
@@ -3441,100 +3180,6 @@
       END SUBROUTINE spinsq
 
 !
-!-------------------------------------------------------------------------
-      SUBROUTINE updatc(ccc,x0,phi,bephi,becp,bec,cp)
-!-----------------------------------------------------------------------
-!     input ccc : dt**2/emass (unchanged in output)
-!     input x0  : converged lambdas from ortho-loop (unchanged in output)
-!     input cp  : non-orthonormal cp=c0+dh/dc*ccc
-!     input bec : <cp|beta_i>
-!     input phi 
-!     output cp : orthonormal cp=cp+lambda*phi
-!     output bec: bec=becp+lambda*bephi
-!
-      USE ions_base, ONLY: nsp, na
-      USE io_global, ONLY: stdout
-      USE cvan, ONLY: nvb, ish
-      USE uspp, ONLY: nhsa => nkb, nhsavb=>nkbus
-      USE uspp_param, ONLY: nh
-      USE gvecw, ONLY: ngw
-      USE electrons_base, ONLY: nx => nbspx, n => nbsp
-      USE control_flags, ONLY: iprint, iprsta
-!
-      IMPLICIT NONE
-!
-      COMPLEX(8) cp(ngw,n), phi(ngw,n)
-      REAL(8)   bec(nhsa,n), x0(nx,nx), ccc
-      REAL(8)   bephi(nhsa,n), becp(nhsa,n)
-! local variables
-      INTEGER i, j, ig, is, iv, ia, inl
-      REAL(8) wtemp(n,nhsa) ! automatic array
-      COMPLEX(8), ALLOCATABLE :: wrk2(:,:)
-!
-!     lagrange multipliers
-!
-      CALL start_clock( 'updatc' )
-      
-      ALLOCATE( wrk2( ngw, n ) )
-
-      wrk2 = (0.d0, 0.d0)
-      DO j=1,n
-         CALL DSCAL(n,ccc,x0(1,j),1)
-      END DO
-!
-!     wrk2 = sum_m lambda_nm s(r(t+dt))|m>
-!
-      CALL MXMA(phi,1,2*ngw,x0,nx,1,wrk2,1,2*ngw,2*ngw,n,n)
-!
-      DO i=1,n
-         DO ig=1,ngw
-            cp(ig,i)=cp(ig,i)+wrk2(ig,i)
-         END DO
-      END DO
-!    
-!     updating of the <beta|c(n,g)>
-!
-!     bec of vanderbilt species are updated 
-!
-      IF(nvb.GT.0)THEN
-         CALL MXMA(x0,1,nx,bephi,nhsa,1,wtemp,1,n,n,n,nhsavb)
-!
-         DO i=1,n
-            DO inl=1,nhsavb
-               bec(inl,i)=wtemp(i,inl)+becp(inl,i)
-            END DO
-         END DO
-      ENDIF
-!
-      IF (iprsta.GT.2) THEN
-         WRITE( stdout,*)
-         DO is=1,nsp
-            IF(nsp.GT.1) THEN
-               WRITE( stdout,'(33x,a,i4)') ' updatc: bec (is)',is
-               WRITE( stdout,'(8f9.4)')                                       &
-     &            ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)
-            ELSE
-               DO ia=1,na(is)
-                  WRITE( stdout,'(33x,a,i4)') ' updatc: bec (ia)',ia
-                  WRITE( stdout,'(8f9.4)')                                    &
-     &            ((bec(ish(is)+(iv-1)*na(is)+ia,i),iv=1,nh(is)),i=1,n)
-               END DO
-            END IF
-            WRITE( stdout,*)
-         END DO
-      ENDIF
-!
-      DO j=1,n
-         CALL DSCAL(n,1.0/ccc,x0(1,j),1)
-      END DO
-
-      DEALLOCATE( wrk2 )
-!
-      CALL stop_clock( 'updatc' )
-!
-      RETURN
-      END SUBROUTINE updatc
-!
 !-----------------------------------------------------------------------
       SUBROUTINE vofrho(nfi,rhor,rhog,rhos,rhoc,tfirst,tlast,           &
      &     ei1,ei2,ei3,irb,eigrb,sfac,tau0,fion)
@@ -3551,7 +3196,7 @@
 !     rhos output: total potential on smooth real space grid
 !
       USE kinds, ONLY: dp
-      USE control_flags, ONLY: iprint, tvlocw, iprsta, thdyn, tpre, tfor, tprnfor
+      USE control_flags, ONLY: iprint, iprsta, thdyn, tpre, tfor, tprnfor
       USE io_global, ONLY: stdout
       USE parameters, ONLY: natx, nsx
       USE ions_base, ONLY: nas => nax, nsp, na, nat
@@ -3841,14 +3486,6 @@
 !
       etot=ekin+eht+epseu+enl+exc+ebac
       IF(tpre) detot=dekin+dh+dps+denl+dxc+dsr
-!
-      IF(tvlocw.AND.tlast)THEN
-#ifdef __PARA
-         CALL write_rho(46,nspin,rhor)
-#else
-         WRITE(46) ((rhor(ir,iss),ir=1,nnr),iss=1,nspin)
-#endif
-      ENDIF
 !
       DEALLOCATE(rhotmp)
       DEALLOCATE(vtemp)

CPV/cplib_meta.f90

@@ -75,7 +75,7 @@
       use cell_base
       use smooth_grid_dimensions, only: nr1s, nr2s, nr3s, &
             nr1sx, nr2sx, nr3sx, nnrsx
-      use electrons_base, only: nx => nbspx, n => nbsp, f, ispin => fspin, nspin
+      use electrons_base, only: nx => nbspx, n => nbsp, f, ispin, nspin
       use constants, only: pi, fpi
 !
       use cdvan

CPV/cpr.f90

@@ -28,7 +28,7 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
   USE energies,                 ONLY : eht, epseu, exc, etot, eself, enl, &
                                        ekin, atot, entropy, egrand, enthal, &
                                        ekincm, print_energies
-  USE electrons_base,           ONLY : nbspx, nbsp, ispin => fspin, f, nspin
+  USE electrons_base,           ONLY : nbspx, nbsp, ispin, f, nspin
   USE electrons_base,           ONLY : nel, iupdwn, nupdwn, nudx, nelt
   USE efield_module,            ONLY : efield, epol, tefield, allocate_efield, &
                                        efield_update, ipolp, qmat, gqq,        &
@@ -130,6 +130,8 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
   USE metadyn_base,             ONLY : set_target
   USE autopilot,                ONLY : pilot
   USE ions_nose,                ONLY : ions_nose_allocate, ions_nose_shiftvar
+  USE orthogonalize,            ONLY : ortho
+  USE orthogonalize_base,       ONLY : updatc
   !
   IMPLICIT NONE
   !
@@ -396,8 +398,8 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
         !
         IF ( tortho ) THEN
            !
-           CALL ortho( eigr, cm, phi, lambda, bigr, iter, ccc, &
+           CALL ortho( eigr, cm(:,:,1,1), phi(:,:,1,1), lambda, bigr, iter, ccc, &
-                       ortho_eps, ortho_max, delt, bephi, becp )
+                       bephi, becp )
            !
         ELSE
            !
@@ -417,7 +419,9 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
         !
         IF ( iprsta >= 3 ) CALL print_lambda( lambda, nbsp, 9, 1.D0 )
         !
-        IF ( tortho ) CALL updatc( ccc, lambda, phi, bephi, becp, bec, cm )
+        IF ( tortho ) &
+           CALL updatc( ccc, nbsp, lambda, SIZE(lambda,1), phi, SIZE(phi,1), &
+                        bephi, SIZE(bephi,1), becp, bec, cm )
         !
         CALL calbec( nvb+1, nsp, eigr, cm, bec )
         !

CPV/electrons.f90

@@ -687,7 +687,7 @@
    SUBROUTINE cp_eigs( nfi, bec, c0, irb, eigrb, rhor, rhog, rhos, lambdap, lambda, tau0, h )
 
      use ensemble_dft, only: tens, ismear, z0, c0diag, becdiag
-     use electrons_base, only: nx => nbspx, n => nbsp, ispin => fspin, f, nspin
+     use electrons_base, only: nx => nbspx, n => nbsp, ispin, f, nspin
      use electrons_base, only: nel, iupdwn, nupdwn, nudx, nelt
      use energies, only: enl, ekin
      use uspp, only: rhovan => becsum

CPV/emptystates.f90

@@ -455,7 +455,7 @@
       COMPLEX(DP), INTENT(INOUT) ::  c_occ(:,:,:,:), c_emp(:,:,:,:)
       TYPE (wave_descriptor), INTENT(IN) ::  wfill, wempt
       TYPE (atoms_type), INTENT(INOUT) :: atoms ! ions structure
-      REAL (DP), INTENT(IN) ::   vpot(:,:,:,:)
+      REAL (DP), INTENT(IN) ::   vpot(:,:)
       LOGICAL, INTENT(IN) :: tortho
       COMPLEX(DP) :: eigr(:,:)
 !
@@ -525,7 +525,7 @@
             CALL nlsm1 ( n_emp, 1, nspnl, eigr, c_emp( 1, 1, ik, ispin ), bece( 1, (ispin-1)*n_emp + 1 ) )
 
             CALL dforce_all( ispin, c_emp(:,:,1,ispin), wempt, fi(:,1,ispin), eforce(:,:,1,ispin), &
-              vpot(:,:,:,ispin), eigr, bece )
+              vpot(:,ispin), eigr, bece )
 
             ! ...       Steepest descent
             DO i = 1, n_emp
@@ -613,7 +613,7 @@
 
       COMPLEX(DP), INTENT(inout) ::  c_emp(:,:,:,:)
       TYPE (wave_descriptor), INTENT(IN) ::  wempt
-      REAL (DP), INTENT(in) ::  vpot(:,:,:,:), fi(:,:,:)
+      REAL (DP), INTENT(in) ::  vpot(:,:), fi(:,:,:)
       COMPLEX (DP) ::  eforce(:,:,:,:)
       LOGICAL, INTENT(IN) :: TORTHO
       COMPLEX(DP) :: eigr(:,:)
@@ -649,7 +649,7 @@
         ! ...   Calculate | dH / dpsi(j) >
         !
         CALL dforce_all( ispin, c_emp(:,:,1,ispin), wempt, fi(:,1,ispin), eforce(:,:,1,ispin), &
-          vpot(:,:,:,ispin), eigr, bece )
+          vpot(:,ispin), eigr, bece )
 
         ! ...     Calculate Eij = < psi(i) | H | psi(j) > = < psi(i) | dH / dpsi(j) >
         DO i = 1, n_emp

CPV/exch_corr.f90

@@ -30,37 +30,28 @@
         SUBROUTINE v2gc( v2xc, grho, rhoer, vpot )
 
           USE kinds,              ONLY: DP
-          USE fft
+          USE fft,                ONLY: pfwfft, pinvfft
-          USE fft_base, ONLY: dfftp
           USE cell_base,          ONLY: tpiba
-          USE mp_global
           USE reciprocal_vectors, ONLY: gstart, gx
+          use grid_dimensions,    only: nnrx
           USE gvecp,              ONLY: ngm
 !                                                                       
           implicit none
 !                                                                       
-          REAL(DP) ::  vpot(:,:,:,:)
+          REAL(DP) ::  vpot(:,:)
-          REAL(DP), intent(in)  ::  v2xc(:,:,:,:,:)
+          REAL(DP), intent(in)  ::  v2xc(:,:,:)
-          REAL(DP), intent(in)  ::  grho(:,:,:,:,:)
+          REAL(DP), intent(in)  ::  grho(:,:,:)
-          REAL(DP), intent(in)  ::  rhoer(:,:,:,:)
+          REAL(DP), intent(in)  ::  rhoer(:,:)
 !                                                                       
-          integer :: ig, ipol, nxl, nyl, nzl, i, j, k, is, js, nspin
+          integer :: ig, ipol, is, js, nspin
-          integer :: ldx, ldy, ldz
+          COMPLEX(DP), allocatable ::  psi(:)
-          COMPLEX(DP), allocatable ::  psi(:,:,:)
           COMPLEX(DP), allocatable ::  vtemp(:)
           COMPLEX(DP), allocatable ::  vtemp_pol(:)
-          REAL(DP), ALLOCATABLE :: v(:,:,:)
+          REAL(DP), ALLOCATABLE :: v(:)
           REAL(DP) :: fac
 ! ...                                                                   
-          ldx   = dfftp%nr1x
+          nspin = SIZE(rhoer,2)
-          ldy   = dfftp%nr2x
-          ldz   = dfftp%npl
-          nxl   = MIN( dfftp%nr1, SIZE( grho, 1 ) )
-          nyl   = MIN( dfftp%nr2, SIZE( grho, 2 ) )
-          nzl   = MIN( dfftp%npl, SIZE( grho, 3 ) )
-          nspin = SIZE(rhoer,4)
 
-          !fac = REAL(nspin)
           fac = 1.0d0
 
           ALLOCATE( vtemp( ngm ) )
@@ -68,20 +59,14 @@
 
           DO js = 1, nspin
             !
-            ALLOCATE( psi( ldx, ldy, ldz ) )
+            ALLOCATE( psi( nnrx ) )
             !
             vtemp = 0.0d0
 
             DO ipol = 1, 3
               DO is = 1, nspin
                 !
-                DO k = 1, nzl
+                psi( 1:nnrx ) = fac * v2xc( 1:nnrx, js, is ) * grho( 1:nnrx, ipol, is )
-                  DO j = 1, nyl
-                    DO i = 1, nxl
-                      psi(i,j,k) = fac * v2xc(i,j,k,js,is) * grho(i,j,k,ipol,is)
-                    END DO
-                  END DO
-                END DO
                 !
                 CALL pfwfft( vtemp_pol, psi )
                 !
@@ -94,17 +79,12 @@
             !
             DEALLOCATE( psi )
 
-            ALLOCATE( v( ldx, ldy, ldz ) )
+            ALLOCATE( v( nnrx ) )
+            v( 1:nnrx ) = 0.0d0
             !
             CALL pinvfft( v, vtemp )
 
-            DO k = 1, nzl
+            vpot( 1:nnrx, js ) = vpot( 1:nnrx, js) - v( 1:nnrx )
-              DO j = 1, nyl
-                DO i = 1, nxl
-                  vpot(i,j,k,js) = vpot(i,j,k,js) - v(i,j,k)
-                END DO
-              END DO
-            END DO
 
             DEALLOCATE( v )
 
@@ -120,25 +100,21 @@
 
     SUBROUTINE stress_gc(grho, v2xc, gcpail, omega)
 !
-      use grid_dimensions, only: nr1, nr2, nr3
+      use grid_dimensions, only: nr1, nr2, nr3, nnrx
-      USE fft_base, ONLY: dfftp
 
         IMPLICIT NONE
 !
-        REAL(DP) ::  v2xc(:,:,:,:,:)
+        REAL(DP) ::  v2xc(:,:,:)
-        REAL(DP) ::  grho(:,:,:,:,:)
+        REAL(DP) ::  grho(:,:,:)
         REAL(DP) ::  gcpail(6)
         REAL(DP) ::  omega
 !
         REAL(DP) :: stre, grhoi, grhoj
-        INTEGER :: i, j, k, ipol, jpol, ic, nxl, nyl, nzl, is, js, nspin
+        INTEGER :: i, ipol, jpol, ic, is, js, nspin
         INTEGER, DIMENSION(6), PARAMETER :: alpha = (/ 1,2,3,2,3,3 /)
         INTEGER, DIMENSION(6), PARAMETER :: beta  = (/ 1,1,1,2,2,3 /)
 ! ...
-        nxl   = MIN( dfftp%nr1, SIZE( grho, 1 ) )
+        nspin = SIZE(grho,3)
-        nyl   = MIN( dfftp%nr2, SIZE( grho, 2 ) )
-        nzl   = MIN( dfftp%npl, SIZE( grho, 3 ) )
-        nspin = SIZE(grho,5)
 
         DO ic = 1, 6
           ipol = alpha(ic)
@@ -146,12 +122,8 @@
           stre = 0.0d0
           DO is = 1, nspin
             DO js = 1, nspin
-              DO k = 1, nzl
+              DO i = 1, nnrx
-                DO j = 1, nyl
+                stre = stre + v2xc(i,is,js) * grho(i,ipol,js) * grho(i,jpol,is)
-                  DO i = 1, nxl
-                    stre = stre + v2xc(i,j,k,is,js) * grho(i,j,k,ipol,js) * grho(i,j,k,jpol,is)
-                  END DO
-                END DO
               END DO
             END DO
           END DO
@@ -184,8 +156,8 @@
       COMPLEX(DP) :: vxc(:,:)
       COMPLEX(DP), INTENT(IN) :: sfac(:,:)
       REAL(DP) :: dexc(:), strvxc
-      REAL(DP) :: grho(:,:,:,:,:)
+      REAL(DP) :: grho(:,:,:)
-      REAL(DP) :: v2xc(:,:,:,:,:)
+      REAL(DP) :: v2xc(:,:,:)
       REAL(DP) :: GAgx_L(:,:)
       REAL(DP) :: rhocp(:,:)
 
@@ -276,51 +248,30 @@
      SUBROUTINE exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxc, vxc, v2xc)
 
         USE kinds,           ONLY: DP
-        USE grid_dimensions, ONLY: nr1l, nr2l, nr3l
+        use grid_dimensions, only: nnrx
         USE funct, ONLY: dft_is_gradient
 
-        REAL (DP) :: rhoetr(:,:,:,:)
+        REAL (DP) :: rhoetr(:,:)
         COMPLEX(DP) :: rhoetg(:,:)
-        REAL (DP) :: grho(:,:,:,:,:)
+        REAL (DP) :: grho(:,:,:)
-        REAL (DP) :: vpot(:,:,:,:)
+        REAL (DP) :: vpot(:,:)
         REAL (DP) :: sxc              ! E_xc   energy
         REAL (DP) :: vxc              ! SUM ( v(r) * rho(r) )
-        REAL (DP) :: v2xc(:,:,:,:,:)
+        REAL (DP) :: v2xc(:,:,:)
-        REAL (DP) :: ddot
+        !
-
+        REAL (DP), EXTERNAL :: ddot
-        INTEGER :: nspin, nnr, ispin, j, k, i
 
+        INTEGER :: nspin, ispin
         logical :: is_gradient
 
         is_gradient =  dft_is_gradient()
 
+
         !  vpot = vxc(rhoetr); vpot(r) <-- u(r)
 
-          nnr   = SIZE( rhoetr, 1 ) * SIZE( rhoetr, 2 ) * SIZE( rhoetr, 3 )
+        nspin = SIZE( rhoetr, 2 )
-          nspin = SIZE( rhoetr, 4 )
-
         !
-          IF( nnr /= nr3l * nr2l * nr1l ) THEN
+        CALL exch_corr_wrapper( nnrx, nspin, grho(1,1,1), rhoetr(1,1), sxc, vpot(1,1), v2xc(1,1,1) )
-            DO ispin = 1, nspin
-              DO k = 1, SIZE( rhoetr, 3 )
-                DO j = 1, SIZE( rhoetr, 2 )
-                  DO i = 1, SIZE( rhoetr, 1 )
-                    IF( i > nr1l .OR. j > nr2l .OR. k > nr3l ) THEN
-                      rhoetr( i, j, k,    ispin ) = 0.0d0
-                      IF( is_gradient ) THEN
-                        grho  ( i, j, k, :, ispin ) = 0.0d0
-                      END IF
-                    END IF
-                  END DO
-                END DO
-              END DO
-            END DO
-          END IF
-
-          !
-          CALL exch_corr_wrapper( nnr, nspin, grho(1,1,1,1,1), rhoetr(1,1,1,1), &
-                                  sxc, vpot(1,1,1,1), v2xc(1,1,1,1,1) )
-
         !
         IF( dft_is_gradient() ) THEN
           ! ... vpot additional term for gradient correction
@@ -332,12 +283,7 @@
         !
         vxc = 0.0d0
         DO ispin = 1, nspin
-            DO k = 1, nr3l
+           vxc = vxc + DDOT ( nnrx, vpot(1,ispin), 1, rhoetr(1,ispin), 1 )
-              DO j = 1, nr2l
-                vxc = vxc + &
-                DDOT ( nr1l, vpot(1,j,k,ispin), 1, rhoetr(1,j,k,ispin), 1 )
-              END DO
-            END DO
         END DO
 
 

CPV/fft.f90

@@ -154,9 +154,9 @@
 
       IMPLICIT NONE
 
-      COMPLEX(DP), INTENT(INOUT) :: cpsi(:,:,:)
+      COMPLEX(DP), INTENT(INOUT) :: cpsi(:)
       COMPLEX(DP), INTENT(OUT) :: C(:)
-      COMPLEX(DP), ALLOCATABLE :: psi(:,:,:)
+      COMPLEX(DP), ALLOCATABLE :: psi(:)
       COMPLEX(DP), ALLOCATABLE :: zwrk(:) 
       REAL(DP) :: t1
       INTEGER :: ierr
@@ -166,20 +166,16 @@
       IF ( first ) &
         CALL errore( ' pfwfft 2 ', ' fft not initialized ', 1 )
 
-      IF ( SIZE( cpsi, 1 ) /= dfftp%nr1x ) THEN
+      IF ( SIZE( cpsi ) /= dfftp%nnr ) THEN
-        WRITE( stdout, * ) 'Values = ', SIZE( cpsi, 1 ), dfftp%nr1x
+        WRITE( stdout, * ) 'Values = ', SIZE( cpsi ), dfftp%nnr
         CALL errore( ' pfwfft 2 ', ' inconsistent array dimensions ', 1 )
       END IF
-      IF ( SIZE( cpsi, 2 ) /= dfftp%nr2x ) &
-        CALL errore( ' pfwfft 2 ', ' inconsistent array dimensions ', 2 )
-      IF ( SIZE( cpsi, 3 ) /= dfftp%npl ) &
-        CALL errore( ' pfwfft 2 ', ' inconsistent array dimensions ', 3 )
 
 #if defined __PARA
 
       ALLOCATE( zwrk( dfftp%nsp( mpime + 1 ) * dfftp%nr3x ) )
 
-      CALL pc3fft_drv(cpsi(1,1,1), zwrk, -1, dfftp, FFT_MODE_POTE)
+      CALL pc3fft_drv(cpsi(1), zwrk, -1, dfftp, FFT_MODE_POTE)
 
       CALL psi2c( zwrk, SIZE( zwrk ), c(1), c(1), ng, 10 )
 
@@ -187,7 +183,7 @@
 
 #else
 
-      ALLOCATE( psi( SIZE( cpsi, 1 ),  SIZE( cpsi, 2 ), SIZE( cpsi, 3 ) ) )
+      ALLOCATE( psi( SIZE( cpsi ) ) )
 
       psi = cpsi
 
@@ -221,10 +217,10 @@
 
       IMPLICIT NONE
 
-      REAL(DP),  INTENT(IN) :: A(:,:,:)
+      REAL(DP),  INTENT(IN) :: A(:)
       COMPLEX(DP) :: C(:)
 
-      COMPLEX(DP), allocatable :: psi(:,:,:)
+      COMPLEX(DP), allocatable :: psi(:)
       COMPLEX(DP), ALLOCATABLE :: zwrk(:) 
       REAL(DP) :: t1
       INTEGER :: ierr, ig, k, is
@@ -234,14 +230,10 @@
       IF ( first ) &
         CALL errore( ' pfwfft 1 ', ' fft not initialized ', 1 )
 
-      IF ( SIZE( A, 1 ) /= dfftp%nr1x ) &
+      IF ( SIZE( A ) /= dfftp%nnr ) &
         CALL errore( ' pfwfft 1 ', ' inconsistent array dimensions ', 1 )
-      IF ( SIZE( A, 2 ) /= dfftp%nr2x ) &
-        CALL errore( ' pfwfft 1 ', ' inconsistent array dimensions ', 2 )
-      IF ( SIZE( A, 3 ) /= dfftp%npl ) &
-        CALL errore( ' pfwfft 1 ', ' inconsistent array dimensions ', 3 )
 
-      ALLOCATE( psi( SIZE(A,1), SIZE(A,2), SIZE(A,3) ), STAT=ierr)
+      ALLOCATE( psi( SIZE(A) ), STAT=ierr)
       IF( ierr /= 0 )  call errore(' PFWFFT ', ' allocation of psi failed ' ,0)
 
       psi = CMPLX( A, 0.d0 )
@@ -250,7 +242,7 @@
 
       ALLOCATE( zwrk( dfftp%nsp( mpime + 1 ) * dfftp%nr3x ) )
 
-      CALL pc3fft_drv(psi(1,1,1), zwrk, -1, dfftp, FFT_MODE_POTE)
+      CALL pc3fft_drv(psi(1), zwrk, -1, dfftp, FFT_MODE_POTE)
 
       CALL psi2c( zwrk(1), SIZE( zwrk ), c(1), c(1), ng, 10 )
 
@@ -287,12 +279,12 @@
 
       IMPLICIT NONE
 
-      REAL(DP), INTENT(INOUT)  :: C(:,:,:)
+      REAL(DP), INTENT(INOUT)  :: C(:)
       REAL(DP), INTENT(IN), OPTIONAL :: ALPHA
       COMPLEX(DP), INTENT(IN) :: A(:)
 
       INTEGER :: ierr
-      COMPLEX(DP), ALLOCATABLE :: psi(:,:,:)
+      COMPLEX(DP), ALLOCATABLE :: psi(:)
       COMPLEX(DP), ALLOCATABLE :: zwrk(:) 
       REAL(DP) t1
 !
@@ -301,14 +293,10 @@
       IF ( first ) &
         CALL errore(' pinvfft ',' fft not initialized ', 0 )
 
-      IF ( SIZE( c, 1 ) /= dfftp%nr1x ) &
+      IF ( SIZE( c ) /= dfftp%nnr ) &
         CALL errore( ' pinvfft 2 ', ' inconsistent array dimensions ', 1 )
-      IF ( SIZE( c, 2 ) /= dfftp%nr2x ) &
-        CALL errore( ' pinvfft 2 ', ' inconsistent array dimensions ', 2 )
-      IF ( SIZE( c, 3 ) /= dfftp%npl ) &
-        CALL errore( ' pinvfft 2 ', ' inconsistent array dimensions ', 3 )
 
-      ALLOCATE( psi( SIZE( c, 1 ), SIZE( c, 2 ), SIZE( c, 3 ) ), STAT=ierr)
+      ALLOCATE( psi( SIZE( c ) ), STAT=ierr)
       IF( ierr /= 0 )  call errore(' PFWFFT ', ' allocation of psi failed ' ,0)
 
 #if defined __PARA
@@ -321,7 +309,7 @@
         CALL c2psi( zwrk, SIZE( zwrk ), a(1), a(1), ng, 11 )
       END IF
 
-      CALL pc3fft_drv(psi(1,1,1), zwrk, +1, dfftp, FFT_MODE_POTE)
+      CALL pc3fft_drv(psi(1), zwrk, +1, dfftp, FFT_MODE_POTE)
 
       DEALLOCATE( zwrk )
 
@@ -370,10 +358,10 @@
 
      COMPLEX(DP) :: C(:)
      COMPLEX(DP), OPTIONAL :: CA(:)
-     COMPLEX(DP) :: psi(:,:,:)
+     COMPLEX(DP) :: psi(:)
      REAL(DP)  :: T1
      INTEGER :: ierr
-     COMPLEX(DP), ALLOCATABLE :: psitmp(:,:,:)
+     COMPLEX(DP), ALLOCATABLE :: psitmp(:)
      COMPLEX(DP), ALLOCATABLE :: zwrk(:) 
 
      t1 = cclock()
@@ -381,19 +369,15 @@
      IF ( first ) &
        CALL errore(' pw_fwfft 1 ',' fft not initialized ', 1 )
 
-     IF ( SIZE( psi, 1 ) /= dffts%nr1x ) &
+     IF ( SIZE( psi ) /= dffts%nnr ) &
        CALL errore( ' pw_fwfft 1 ', ' inconsistent array dimensions ', 1 )
-     IF ( SIZE( psi, 2 ) /= dffts%nr2x ) &
-       CALL errore( ' pw_fwfft 1 ', ' inconsistent array dimensions ', 2 )
-     IF ( SIZE( psi, 3 ) /= dffts%npl ) &
-       CALL errore( ' pw_fwfft 1 ', ' inconsistent array dimensions ', 3 )
 
 #if defined __PARA
 
 
      ALLOCATE( zwrk( dffts%nsp( mpime + 1 ) * dffts%nr3x ) )
 
-     CALL pc3fft_drv(psi(1,1,1), zwrk, -1, dffts, FFT_MODE_WAVE)
+     CALL pc3fft_drv(psi(1), zwrk, -1, dffts, FFT_MODE_WAVE)
 
      IF( PRESENT( ca ) ) THEN
        CALL psi2c( zwrk, SIZE( zwrk ), c(1), ca(1), ngw, 2 )
@@ -405,7 +389,7 @@
 
 #else
 
-     ALLOCATE( psitmp( SIZE( psi, 1 ), SIZE( psi, 2 ), SIZE( psi, 3 ) ) )
+     ALLOCATE( psitmp( SIZE( psi ) ) )
  
      psitmp = psi
 
@@ -441,7 +425,7 @@
 
      COMPLEX(DP), INTENT(IN) :: C(:)
      COMPLEX(DP), INTENT(IN), OPTIONAL :: CA(:)
-     COMPLEX(DP) :: psi(:,:,:)
+     COMPLEX(DP) :: psi(:)
      COMPLEX(DP), ALLOCATABLE :: zwrk(:) 
      REAL(DP) :: T1
 
@@ -451,12 +435,8 @@
      T1 = cclock()
 
 
-     IF ( SIZE( psi, 1 ) /= dffts%nr1x ) &
+     IF ( SIZE( psi ) /= dffts%nnr ) &
        CALL errore( ' pw_invfft 1 ', ' inconsistent array dimensions ', 1 )
-     IF ( SIZE( psi, 2 ) /= dffts%nr2x ) &
-       CALL errore( ' pw_invfft 1 ', ' inconsistent array dimensions ', 2 )
-     IF ( SIZE( psi, 3 ) /= dffts%npl ) &
-       CALL errore( ' pw_invfft 1 ', ' inconsistent array dimensions ', 3 )
 
 #if defined __PARA
 
@@ -472,7 +452,7 @@
        END IF
      END IF
 
-     CALL pc3fft_drv(psi(1,1,1), zwrk, +1, dffts, FFT_MODE_WAVE)
+     CALL pc3fft_drv(psi(1), zwrk, +1, dffts, FFT_MODE_WAVE)
 
      DEALLOCATE( zwrk )
 

CPV/fields_type.f90

@@ -9,7 +9,6 @@
 
         USE kinds, ONLY: DP
         USE parallel_types, ONLY: descriptor, processors_grid
-        USE descriptors_module, ONLY: desc_init
         IMPLICIT NONE
         PRIVATE
         SAVE

CPV/forces.f90

@@ -43,13 +43,13 @@
       COMPLEX(DP), INTENT(OUT) :: dco(:), dce(:)
       COMPLEX(DP), INTENT(IN)  :: co(:), ce(:)
       REAL(DP),    INTENT(IN)  :: fio, fie
-      REAL(DP),    INTENT(IN)  :: v(:,:,:)
+      REAL(DP),    INTENT(IN)  :: v(:)
       REAL(DP),    INTENT(IN)  :: hg(:)
-      COMPLEX(DP), OPTIONAL    :: psi_stored(:,:,:)
+      COMPLEX(DP), OPTIONAL    :: psi_stored(:)
 
       ! ... declare other variables
       !
-      COMPLEX(DP), ALLOCATABLE :: psi(:,:,:)
+      COMPLEX(DP), ALLOCATABLE :: psi(:)
       COMPLEX(DP) :: fp, fm, aro, are
       REAL(DP)    :: fioby2, fieby2, arg
       INTEGER      :: ig
@@ -60,7 +60,7 @@
         psi_stored = psi_stored * CMPLX(v, 0.0d0)
         CALL pw_fwfft(psi_stored, dco, dce)
       ELSE
-        ALLOCATE( psi(SIZE(v,1), SIZE(v,2), SIZE(v,3)) )
+        ALLOCATE( psi( SIZE(v) ) )
         CALL pw_invfft(psi, co, ce)
         psi = psi * CMPLX(v, 0.0d0)
         CALL pw_fwfft(psi, dco, dce)
@@ -268,7 +268,7 @@
        INTEGER,      INTENT(IN)  :: ib, iss     ! band and spin index
        COMPLEX(DP), INTENT(IN)  :: c(:,:)
        COMPLEX(DP), INTENT(OUT) :: df(:), da(:)
-       REAL (DP),   INTENT(IN)  :: v(:,:,:), bec(:,:), f(:)
+       REAL (DP),   INTENT(IN)  :: v(:), bec(:,:), f(:)
        COMPLEX(DP), INTENT(IN)  :: eigr(:,:)
        type (wave_descriptor), INTENT(IN) :: cdesc
        !
@@ -327,7 +327,7 @@
         INTEGER,                INTENT(IN)    :: ispin
         COMPLEX(DP),           INTENT(INOUT) :: c(:,:)
         type (wave_descriptor), INTENT(IN)    :: cdesc
-        REAL(DP),              INTENT(IN)    :: vpot(:,:,:), f(:)
+        REAL(DP),              INTENT(IN)    :: vpot(:), f(:)
         COMPLEX(DP),           INTENT(OUT)   :: cgrad(:,:)
         COMPLEX(DP),           INTENT(IN)    :: eigr(:,:)
         REAL(DP),              INTENT(IN)    :: bec(:,:)

CPV/fromscra.f90

@@ -25,7 +25,7 @@ MODULE from_scratch_module
   CONTAINS
   !
   !--------------------------------------------------------------------------
-  SUBROUTINE from_scratch_fpmd( rhoe, desc, cm, c0, cp, ce, cdesc, edesc, &
+  SUBROUTINE from_scratch_fpmd( rhoe, cm, c0, cp, ce, cdesc, edesc, &
                                 eigr, ei1, ei2, ei3, sfac, fi, ht, atoms, &
                                 bec, becdr, vpot, edft )
     !------------------------------------------------------------------------
@@ -47,7 +47,6 @@ MODULE from_scratch_module
     USE orthogonalize,    ONLY : ortho
     USE control_flags,    ONLY : tcarpar, tfor, thdyn, tortho, tpre, tranp, &
                                  force_pairing, iprsta, tprnfor, amprp, tsde
-    USE charge_types,     ONLY : charge_descriptor
     USE time_step,        ONLY : delt
     USE runcp_module,     ONLY : runcp_ncpp
     use grid_dimensions,  only : nr1, nr2, nr3
@@ -69,15 +68,14 @@ MODULE from_scratch_module
     COMPLEX(DP),             INTENT(OUT) :: ei2(:,:)
     COMPLEX(DP),             INTENT(OUT) :: ei3(:,:)
     COMPLEX(DP),             INTENT(OUT) :: sfac(:,:)
-    REAL(DP),                INTENT(OUT) :: rhoe(:,:,:,:)
+    REAL(DP),                INTENT(OUT) :: rhoe(:,:)
     REAL(DP),                INTENT(OUT) :: bec(:,:)
     REAL(DP),                INTENT(OUT) :: becdr(:,:,:)
     REAL(DP),                INTENT(OUT) :: fi(:,:,:)
-    REAL(DP),                INTENT(OUT) :: vpot(:,:,:,:)
+    REAL(DP),                INTENT(OUT) :: vpot(:,:)
     TYPE(atoms_type)    ,    INTENT(OUT) :: atoms
     TYPE(dft_energy_type) ,  INTENT(OUT) :: edft
     TYPE(boxdimensions)  ,   INTENT(INOUT) :: ht 
-    TYPE(charge_descriptor), INTENT(IN)    :: desc
     TYPE(wave_descriptor),   INTENT(IN)    :: cdesc, edesc
     COMPLEX(DP),             INTENT(INOUT) :: cm(:,:,:,:), c0(:,:,:,:)
     COMPLEX(DP),             INTENT(INOUT) :: cp(:,:,:,:), ce(:,:,:,:)
@@ -164,9 +162,9 @@ MODULE from_scratch_module
     !
     edft%enl = nlrh_m( cm, cdesc, ttforce, atoms, fi, bec, becdr, eigr )
     !
-    CALL rhoofr( 0, cm, cdesc, fi, rhoe, desc, ht )
+    CALL rhoofr( 0, cm, cdesc, fi, rhoe, ht )
     !
-    CALL vofrhos( ttprint, ttforce, tstress, rhoe, desc, atoms, &
+    CALL vofrhos( ttprint, ttforce, tstress, rhoe, atoms, &
                   vpot, bec, cm, cdesc, fi, eigr, ei1, ei2, ei3, &
                   sfac, timepre, ht, edft )
     !
@@ -255,6 +253,8 @@ MODULE from_scratch_module
     USE runcp_module,         ONLY : runcp_uspp
     USE electrons_base,       ONLY : f, nspin
     USE phase_factors_module, ONLY : strucf
+    USE orthogonalize,        ONLY : ortho
+    USE orthogonalize_base,   ONLY : updatc, calphi
     !
     IMPLICIT NONE
     !
@@ -394,11 +394,10 @@ MODULE from_scratch_module
 !     calphi calculates phi
 !     the electron mass rises with g**2
 !
-      CALL calphi( cm, ngw, ema0bg, bec, nkb, vkb, phi, nbsp )
+      CALL calphi( cm, ngw, bec, nkb, vkb, phi, nbsp, ema0bg )
-
 
       if( tortho ) then
-         CALL ortho( eigr, c0, phi, lambda, bigr, iter, ccc, ortho_eps, ortho_max, delt, bephi, becp )
+         CALL ortho( eigr, c0(:,:,1,1), phi(:,:,1,1), lambda, bigr, iter, ccc, bephi, becp )
       else
          CALL gram( vkb, bec, nkb, c0, ngw, nbsp )
       endif
@@ -410,8 +409,11 @@ MODULE from_scratch_module
 
       if ( tpre ) CALL nlfh( bec, dbec, lambda )
 !
-      if ( tortho ) CALL updatc( ccc, lambda, phi, bephi, becp, bec, c0 )
+      if ( tortho ) CALL updatc( ccc, nbsp, lambda, SIZE(lambda,1), phi, SIZE(phi,1), &
+                                 bephi, SIZE(bephi,1), becp, bec, c0 )
+      !
       CALL calbec ( nvb+1, nsp, eigr, c0, bec )
+
       if ( tpre ) CALL caldbec( ngw, nkb, nbsp, 1, nsp, eigr, cm, dbec, .true. )
 
       if(iprsta.ge.3) CALL dotcsc(eigr,c0)

CPV/grid.f90

@@ -82,34 +82,26 @@
             RETURN
           END SUBROUTINE free_blacs_grid
 
-          SUBROUTINE get_blacs_grid(grid, rows, columns, debug)
+
+
+          SUBROUTINE get_blacs_grid(grid, debug)
+
             TYPE (processors_grid), INTENT(OUT) :: grid
-            INTEGER, INTENT(IN), OPTIONAL :: rows
-            INTEGER, INTENT(IN), OPTIONAL :: columns
             INTEGER, INTENT(IN), OPTIONAL :: debug
             INTEGER :: iam, nproc , nprow, npcol, context, myrow, mycol 
 
+            INTEGER :: ndims, dims(2), coor(2)
+            LOGICAL :: periods(2), reorder
+            INTEGER :: comm_cart
+            INTEGER :: ierr
+
 #if defined __SCALAPACK
             CALL BLACS_PINFO( iam, nproc  )
 #else
-            nproc  = -1
+            ndims = 2
+
 #endif
 
-            IF(.NOT.PRESENT(rows) .AND. .NOT.PRESENT(columns) ) THEN
-              CALL calculate_grid_dims(nproc , nprow, npcol)
-            ELSE IF (PRESENT(rows) .AND. .NOT.PRESENT(columns) ) THEN
-              !IF( rows .GT. nproc  ) THEN
-              !END IF 
-              nprow = rows; npcol = nproc  / rows
-            ELSE IF (.NOT.PRESENT(rows) .AND. PRESENT(columns) ) THEN
-              !IF( columns .GT. nproc  ) THEN
-              !END IF 
-              npcol = columns; nprow = nproc  / columns
-            ELSE
-              !IF( rows * columns .GT. nproc  ) THEN
-              !END IF 
-              nprow = rows; npcol = columns
-            END IF
 
 #if defined __SCALAPACK
             CALL BLACS_GET( -1, 0, context )

CPV/guess.f90

@@ -6,36 +6,12 @@
 ! or http://www.gnu.org/copyleft/gpl.txt .
 !
 #include "f_defs.h"
-#if defined __T3E
-#  define daxpy saxpy
-#  define zaxpy caxpy
-#endif
-
-!  AB INITIO COSTANT PRESSURE MOLECULAR DYNAMICS
-!  ----------------------------------------------
-!  Car-Parrinello Parallel Program
-!  Carlo Cavazzoni - Gerardo Ballabio
-!  SISSA, Trieste, Italy - 1997-99
-!  Last modified: Sun Nov 21 11:19:43 MET 1999
-!  ----------------------------------------------
-!  BEGIN manual
 
 MODULE guess
 
-!  (describe briefly what this module does...)
-!  ----------------------------------------------
-!  routines in this module:
-!  SUBROUTINE guess_setup(tguess_inp)
-!  SUBROUTINE guessc0(tk,c0,cm)
-!  SUBROUTINE guessrho(rho,cm,c0,occ,ht) 
-!  SUBROUTINE ucalc_kp(a,b,ngw,gzero,n,lambda)
-!  SUBROUTINE ucalc(a,b,ngw,gzero,n,lambda)
-!  ----------------------------------------------
-!  END manual
-
 ! ...   declare modules
         USE kinds
-        USE parallel_toolkit, ONLY: matmulp, cmatmulp, &
+        USE parallel_toolkit, ONLY: rep_matmul_drv, &
           diagonalize, cdiagonalize
 
         IMPLICIT NONE
@@ -43,7 +19,7 @@
 
         PRIVATE
 
-        REAL(DP), ALLOCATABLE :: rho_save( :, :, :, : )
+        REAL(DP), ALLOCATABLE :: rho_save( :, : )
    
 ! ...   declare module-scope variables
         LOGICAL :: tguess
@@ -140,7 +116,8 @@
               DO ik = 1, nk
 
                 CALL ucalc_kp(cm(:,:,ik,1),c0(:,:,ik,1),ngw,cdesc%gzero,n,cuu)
-                CALL cmatmulp('N','C',cuu,cuu,ca,n)
+                ! CALL cmatmulp('N','C',cuu,cuu,ca,n)
+                CALL errore(' guess ', ' complex matrix mult to be implemented ', 1 )
                 CALL cdiagonalize(1,ca,costemp,cap,n,nproc,mpime)
                 DO j=1,n
                   DO i=1,n
@@ -150,7 +127,8 @@
                 DO i=1,n
                   costh2(i)=1.0d0/sqrt(costemp(n-i+1))
                 END DO
-                CALL cmatmulp('N','N',cuu,ca,cap,n)
+                !CALL cmatmulp('N','N',cuu,ca,cap,n)
+                CALL errore(' guess ', ' complex matrix mult to be implemented ', 1 )
                 DO j=1,n
                   DO i=1,n
                     cap(i,j)=cap(i,j) * costh2(i)
@@ -186,7 +164,7 @@
               ALLOCATE(crot(ngw,n))
 
               CALL ucalc(cm(:,:,1,1),c0(:,:,1,1),ngw,cdesc%gzero,n,uu)
-              CALL matmulp('T','N',uu,uu,a,n)
+              CALL rep_matmul_drv('T','N',n,n,n,1.0d0,uu,n,uu,n,0.0d0,a,n,group)
               CALL diagonalize(1,a,costemp,ap,n,nproc,mpime)
               DO j=1,n
                 DO i=1,n
@@ -196,7 +174,7 @@
               DO i=1,n
                 costh2(i)=1.0d0/sqrt(costemp(n-i+1))
               END DO
-              CALL matmulp('N','N',uu,a,ap,n)
+              CALL rep_matmul_drv('N','N',n,n,n,1.0d0,uu,n,a,n,0.0d0,ap,n,group)
               DO j=1,n
                 DO i=1,n
                   ap(i,j)=ap(i,j) * costh2(i)
@@ -245,7 +223,7 @@
       
 !  ----------------------------------------------
 !  ----------------------------------------------
-          SUBROUTINE guessrho(rho, desc, cm, c0, cdesc, occ, ht ) 
+          SUBROUTINE guessrho(rho, cm, c0, cdesc, occ, ht ) 
 
 !  (describe briefly what this routine does...)
 !  ----------------------------------------------
@@ -255,38 +233,33 @@
              use brillouin, only: kpoints, kp
              USE wave_types
              USE parameters, ONLY: nspinx
-             USE charge_types, ONLY: charge_descriptor
 
 ! ...        declare subroutine argument
-             REAL(DP), INTENT(OUT) :: rho(:,:,:,:)
+             REAL(DP), INTENT(OUT) :: rho(:,:)
-             TYPE (charge_descriptor), INTENT(IN) :: desc
              COMPLEX(DP), INTENT(IN) :: c0(:,:,:,:), cm(:,:,:,:)
              TYPE (wave_descriptor), INTENT(IN) :: cdesc
              TYPE (boxdimensions), INTENT(IN) :: ht
              REAL(DP), INTENT(IN) :: occ(:,:,:)
 
 ! ...        declare other variables
-             REAL(DP), ALLOCATABLE :: rho0( :, :, :, : )
+             REAL(DP), ALLOCATABLE :: rho0( :, : )
 
              LOGICAL, SAVE :: tfirst = .TRUE.
-             INTEGER :: ispin, nspin, nx, ny, nz
+             INTEGER :: ispin, nspin
 
 ! ...        end of declarations
 !  ----------------------------------------------
 
-             nx = SIZE( rho, 1 )
+             nspin = SIZE( rho, 2 )
-             ny = SIZE( rho, 2 )
-             nz = SIZE( rho, 3 )
-             nspin = SIZE( rho, 4 )
 
              IF( tfirst ) THEN
-               ALLOCATE( rho_save( nx, ny, nz, nspin ) )
+               ALLOCATE( rho_save( SIZE( rho, 1 ), nspin ) )
-               CALL rhoofr( 1, cm, cdesc, occ, rho_save, desc, ht)
+               CALL rhoofr( 1, cm, cdesc, occ, rho_save, ht)
                tfirst = .FALSE.
              END IF
 
-             ALLOCATE( rho0( nx, ny, nz, nspin ) )
+             ALLOCATE( rho0( SIZE( rho, 1 ), nspin ) )
-             CALL rhoofr( 1, c0, cdesc, occ, rho0, desc, ht)
+             CALL rhoofr( 1, c0, cdesc, occ, rho0, ht)
 
              rho = 2.0d0 * rho0 - rho_save
 

CPV/init_run.f90

@@ -51,8 +51,8 @@ SUBROUTINE init_run()
                                        becdr, sfac, eigr, ei1, ei2, ei3, taub, &
                                        irb, eigrb, rhog, rhos, rhor, bephi,    &
                                        becp, acc, acc_this_run, wfill, wempt,  &
-                                       edft, nfi, atoms0, vpot, occn, desc,    &
+                                       edft, nfi, atoms0, vpot, occn,     &
-                                       rhoe, atomsm, ht0, htm
+                                       atomsm, ht0, htm
   USE cp_main_variables,        ONLY : allocate_mainvar
   USE energies,                 ONLY : eself, enl, ekin, etot, enthal, ekincm
   USE stre,                     ONLY : stress
@@ -257,7 +257,7 @@ SUBROUTINE init_run()
         !
      ELSE IF ( program_name == 'FPMD' ) THEN
         !
-        CALL from_scratch( rhoe, desc, cm, c0, cp, ce, wfill, wempt, eigr, &
+        CALL from_scratch( rhor, cm, c0, cp, ce, wfill, wempt, eigr, &
                            ei1, ei2, ei3, sfac, occn, ht0, atoms0, bec, &
                            becdr, vpot, edft )
         !
@@ -280,7 +280,7 @@ SUBROUTINE init_run()
      ELSE IF( program_name == 'FPMD' ) THEN
         !
         CALL readfile( nfi, tps, c0, cm, wfill, occn, atoms0, atomsm, acc,     &
-                       taui, cdmi, htm, ht0, rhoe, desc, vpot)
+                       taui, cdmi, htm, ht0, rhor, vpot)
         !
      END IF
      !
@@ -294,7 +294,7 @@ SUBROUTINE init_run()
         !
      ELSE IF( program_name == 'FPMD' ) THEN
         !
-        CALL from_restart( nfi, acc, rhoe, desc, cm, c0, wfill, eigr, ei1, ei2, &
+        CALL from_restart( nfi, acc, rhor, cm, c0, wfill, eigr, ei1, ei2, &
                            ei3, sfac, occn, htm, ht0, atomsm, atoms0, bec, &
                            becdr, vpot, edft)
         !

CPV/inner_loop.f90

@@ -24,7 +24,7 @@
       use energies, only: eht, epseu, exc, etot, eself, enl, ekin,          &
      &                    atot, entropy, egrand
       use electrons_base, only: f, nspin, nel, iupdwn, nupdwn, nudx, nelt, &
-                                nx => nbspx, n => nbsp, ispin => fspin
+                                nx => nbspx, n => nbsp, ispin 
 
       use ensemble_dft, only: tens, tgrand, ninner, ismear, etemp, ef,       &
      &                tdynz, tdynf, zmass, fmass, fricz, fricf, z0, c0diag,  &

CPV/input.f90

@@ -320,8 +320,8 @@ MODULE input
           ! gvectors and charge density, in reciprocal space.
      !
      trhor_  = ( TRIM( calculation ) == 'nscf' )
-     trhow_  = ( TRIM( disk_io ) == 'high' )
+     trhow_  = ( TRIM( disk_io ) == 'high' )     !  charge density now written to XML file
-     tvlocw_ = .FALSE.
+     tvlocw_ = ( TRIM( disk_io ) == 'high' )     !  warning this is not working
      !
      SELECT CASE( TRIM( verbosity ) )
        CASE( 'minimal' )
@@ -1354,9 +1354,9 @@ MODULE input
 590 FORMAT( 3X,'Electron temperature control via nose thermostat')
     !
 700 FORMAT( /,3X, 'Verbosity: iprsta = ',i2,/)
-720 FORMAT(   3X, 'charge density is read from unit 47')
+720 FORMAT(   3X, 'charge density is read from file')
-721 FORMAT(   3X, 'charge density is written in unit 47')
+721 FORMAT(   3X, 'warning trhow has no effect rho is now written to XML save file')
-722 FORMAT(   3X, 'local potential is written in unit 46')
+722 FORMAT(   3X, 'warning tvlocw has no effect vloc is not written to file')
     !
   END SUBROUTINE modules_info
   !

CPV/ksstates.f90

@@ -238,7 +238,7 @@
         TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
         COMPLEX(DP)  ::  eigr(:,:)
         REAL(DP), INTENT(IN)  ::  occ(:,:,:), bec(:,:)
-        REAL (DP) ::  vpot(:,:,:,:)
+        REAL (DP) ::  vpot(:,:)
 
 ! ...   declare other variables
         INTEGER ::  i, ik, ib, nk, ig, ngw, nb_g, nb_l, ispin, nspin, iks
@@ -279,7 +279,7 @@
             ALLOCATE(  eforce( ngw,  nb_l, nk ))
 
             CALL dforce_all( ispin, cf(:,:,1,ispin_wfc), wfill, occ(:,1,ispin), eforce(:,:,1), &
-              vpot(:,:,:,ispin), eigr, bec )
+              vpot(:,ispin), eigr, bec )
 
             CALL kohn_sham( ispin, cf(:,:,:,ispin_wfc), wfill, eforce )
 
@@ -302,7 +302,7 @@
               ALLOCATE(  eforce( ngw,  nb_l, nk ))
 
               CALL dforce_all( ispin, ce(:,:,1,ispin), wempt, fi(:,1), eforce(:,:,1), &
-                               vpot(:,:,:,ispin), eigr, bec )
+                               vpot(:,ispin), eigr, bec )
 
               CALL kohn_sham( ispin, ce(:,:,:,ispin), wempt, eforce )
 
@@ -400,7 +400,7 @@
         TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
         COMPLEX(DP)  ::  eigr(:,:)
         REAL(DP), INTENT(IN)  ::  occ(:,:,:), bec(:,:)
-        REAL (DP) ::  vpot(:,:,:,:)
+        REAL (DP) ::  vpot(:,:)
 
 ! ...   declare other variables
         INTEGER ::  i, ik, ib, nk, ig, ngw, nb_g, nb_l, iks, nb, ispin
@@ -440,9 +440,9 @@
           ALLOCATE(  eforce( ngw,  nb, 1, 2 ) )
 
           CALL dforce_all( 1, cf(:,:,1,1), wfill, occ(:,1,1), eforce(:,:,1,1), &
-              vpot(:,:,:,1), eigr, bec )
+              vpot(:,1), eigr, bec )
           CALL dforce_all( 2, cf(:,:,1,1), wfill, occ(:,1,2), eforce(:,:,1,2), &
-              vpot(:,:,:,2), eigr, bec )
+              vpot(:,2), eigr, bec )
 
           DO i = 1, nupdwn(2)
             eforce(:,i,1,1) = occ(i,1,1) * eforce(:,i,1,1) + occ(i,1,2) * eforce(:,i,1,2)
@@ -471,12 +471,12 @@
 
             ALLOCATE(  eforce( ngw,  nb_l, 1, 1 ))
 
-            CALL dforce_all( 1, ce(:,:,1,1), wempt, fi(:,1), eforce(:,:,1,1), vpot(:,:,:,1), &
+            CALL dforce_all( 1, ce(:,:,1,1), wempt, fi(:,1), eforce(:,:,1,1), vpot(:,1), &
                              eigr, bec )
 
             CALL kohn_sham( 1, ce(:,:,:,1), wempt, eforce(:,:,:,1) )
 
-            CALL dforce_all( 2, ce(:,:,1,2), wempt, fi(:,1), eforce(:,:,1,1), vpot(:,:,:,2), &
+            CALL dforce_all( 2, ce(:,:,1,2), wempt, fi(:,1), eforce(:,:,1,1), vpot(:,2), &
                              eigr, bec )
 
             CALL kohn_sham( 2, ce(:,:,:,2), wempt, eforce(:,:,:,1) )
@@ -555,6 +555,7 @@
         USE io_global, ONLY: ionode, ionode_id
         USE io_global, ONLY: stdout
         USE fft_base, ONLY: dfftp
+        USE grid_dimensions, ONLY: nr1, nr2, nr3, nr1x, nr2x, nr3x, nnrx
 
         IMPLICIT NONE
 
@@ -562,27 +563,13 @@
         CHARACTER(LEN=*) :: file_name
         COMPLEX(DP), ALLOCATABLE :: zcomp(:)
         REAL(DP), ALLOCATABLE :: rcomp2(:)
-        COMPLEX(DP), ALLOCATABLE :: psi2(:,:,:)
+        COMPLEX(DP), ALLOCATABLE :: psi2(:)
-        INTEGER   ::  nr1_l, nr2_l, nr3_l, nr1_g, nr2_g, nr3_g
         INTEGER   ::  i, j, k, istr, izl
         REAL(DP) :: charge
         LOGICAL   :: top
 
-        nr1_g = dfftp%nr1
+        ALLOCATE( zcomp( nr3 ), rcomp2( nr3 ) )
-        nr2_g = dfftp%nr2
+        ALLOCATE( psi2( nnrx ) )
-        nr3_g = dfftp%nr3
-
-        nr1_l = dfftp%nr1x
-        nr2_l = dfftp%nr2x
-        nr3_l = dfftp%npl
-
-        izl = 1
-        DO i = 1, mpime
-          izl = izl + dfftp%npp( i )
-        END DO
-
-        ALLOCATE( zcomp( nr3_g ), rcomp2( nr3_g ) )
-        ALLOCATE( psi2( nr1_l, nr2_l, nr3_l ) )
 
         CALL pw_invfft( psi2, psi, psi )
         
@@ -597,25 +584,39 @@
         END IF
 
         charge = 0.0d0
-        DO i = 1, nr1_g
+
-          DO j = 1, nr2_g
+        izl = dfftp%ipp( mpime + 1 )
+
+        DO i = 1, nr1
+
+          DO j = 1, nr2
+
             zcomp = 0.0d0
 
-            zcomp( izl : ( izl + nr3_l - 1 ) ) = psi2( i, j, 1 : nr3_l )
+            istr = i + nr1 * ( j - 1 )
-            CALL mp_sum( zcomp(1:nr3_g) )
+
+            DO k = 1, dfftp%npl 
+               zcomp( izl + k ) = psi2( istr + nr1 * nr2 * ( k - 1 ) )
+            END DO
+
+            CALL mp_sum( zcomp( 1 : nr3 ) )
+
             IF ( ionode ) THEN
               rcomp2 = DBLE(zcomp)**2
-              WRITE(ksunit, fmt='(F10.5)') ( rcomp2(k), k=1, nr3_g )
+              WRITE(ksunit, fmt='(F10.5)') ( rcomp2(k), k=1, nr3 )
               charge = charge + SUM(rcomp2)
             END IF
+
             CALL mp_barrier()
+
           END DO
+
         END DO
 
         IF ( ionode ) THEN
           CLOSE(ksunit)
           WRITE( stdout,'(3X,A15," integrated charge : ",F14.5)')  &
-     &      file_name(1:istr), charge / DBLE(nr1_g*nr2_g*nr3_g)
+     &      file_name(1:istr), charge / DBLE(nr1*nr2*nr3)
         END IF
         DEALLOCATE(zcomp, rcomp2, psi2)
 ! ...

CPV/main.f90

@@ -187,7 +187,7 @@
       USE io_files        , ONLY: outdir, prefix
       USE printout_base   , ONLY: printout_base_init
       USE cp_main_variables, ONLY : atoms0, atomsp, atomsm, ei1, ei2, ei3, eigr, sfac, &
-                                    ht0, htm, htp, rhoe, vpot, desc, wfill, wempt,     &
+                                    ht0, htm, htp, rhor, vpot, wfill, wempt,     &
                                     acc, acc_this_run, occn, edft, nfi, bec, becdr
       USE cg_module,         ONLY : tcg
       IMPLICIT NONE
@@ -344,7 +344,7 @@
            !
            ! ...     perform DIIS minimization on electronic states
            !
-           CALL runsdiis(ttprint, rhoe, desc, atoms0, bec, becdr, &
+           CALL runsdiis(ttprint, rhor, atoms0, bec, becdr, &
                eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht0, occn, ei, &
                vpot, doions, edft )
            !
@@ -354,7 +354,7 @@
            !
            IF( nspin > 1 ) CALL errore(' cpmain ',' lsd+diis not allowed ',0)
            !
-           CALL rundiis(ttprint, rhoe, desc, atoms0, bec, becdr, &
+           CALL rundiis(ttprint, rhor, atoms0, bec, becdr, &
                eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht0, occn, ei, &
                vpot, doions, edft )
            !
@@ -362,7 +362,7 @@
            !
            ! ...     on entry c0 should contain the wavefunctions to be optimized
            !
-           CALL runcg(tortho, ttprint, rhoe, desc, atoms0, bec, becdr, &
+           CALL runcg(tortho, ttprint, rhor, atoms0, bec, becdr, &
                eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht0, occn, ei, &
                vpot, doions, edft, ekin_maxiter, etot_conv_thr, tconv_cg )
            !
@@ -371,13 +371,13 @@
            !
         ELSE IF ( tsteepdesc ) THEN
            !
-           CALL runsd(tortho, ttprint, ttforce, rhoe, desc, atoms0, bec, becdr, &
+           CALL runsd(tortho, ttprint, ttforce, rhor, atoms0, bec, becdr, &
                eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht0, occn, ei, &
                vpot, doions, edft, ekin_maxiter, ekin_conv_thr )
            !
         ELSE IF ( tconjgrad_ion%active ) THEN
            !
-           CALL runcg_ion(nfi, tortho, ttprint, rhoe, desc, atomsp, atoms0, &
+           CALL runcg_ion(nfi, tortho, ttprint, rhor, atomsp, atoms0, &
                atomsm, bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht0, occn, ei, &
                vpot, doions, edft, tconvthrs%derho, tconvthrs%force, tconjgrad_ion%nstepix, &
                tconvthrs%ekin, tconjgrad_ion%nstepex )
@@ -406,9 +406,9 @@
         s5      = cclock()
         timernl = s5 - s4
 
-        ! ...   compute the new charge density "rhoe"
+        ! ...   compute the new charge density "rhor"
         !
-        CALL rhoofr( nfi, c0, wfill, occn, rhoe, desc, ht0)
+        CALL rhoofr( nfi, c0, wfill, occn, rhor, ht0)
 
         IF(memchk) CALL memstat(6)
 
@@ -418,7 +418,7 @@
         ! ...   vofrhos compute the new DFT potential "vpot", and energies "edft",
         ! ...   ionc forces "fion" and stress "pail".
         !
-        CALL vofrhos(ttprint, ttforce, tstress, rhoe, desc, atoms0, &
+        CALL vofrhos(ttprint, ttforce, tstress, rhor, atoms0, &
           vpot, bec, c0, wfill, occn, eigr, ei1, ei2, ei3, sfac, timepre, ht0, edft)
 
         ! CALL debug_energies( edft ) ! DEBUG
@@ -578,7 +578,7 @@
 
         IF( self_interaction /= 0 ) THEN
            IF ( nat_localisation > 0 .AND. ttprint ) THEN
-              CALL localisation( cp( : , nupdwn(1), 1, 1 ), atoms0, ht0, desc)
+              CALL localisation( cp( : , nupdwn(1), 1, 1 ), atoms0, ht0)
            END IF
         END IF
 
@@ -700,7 +700,7 @@
         !
         IF( ttsave .OR. ttexit ) THEN
           CALL writefile( nfi, tps, c0, cm, wfill, occn, atoms0, atomsm, acc,  &
-                          taui, cdmi, htm, ht0, rhoe, desc, vpot )
+                          taui, cdmi, htm, ht0, rhor, vpot )
         END IF
 
         IF( ttexit .AND. .NOT. ttprint ) THEN
@@ -746,12 +746,12 @@
       END IF
 
       IF(tprnsfac) THEN
-        CALL print_sfac(rhoe, desc, sfac)
+        CALL print_sfac(rhor, sfac)
       END IF
 
       ! ... report statistics
 
-      CALL printacc(nfi, rhoe, desc, atomsm, htm, nstep_this_run, acc, acc_this_run)
+      CALL printacc(nfi, nstep_this_run, acc, acc_this_run)
       CALL mp_report_buffers()
       CALL mp_report()
 

CPV/mainvar.f90

@@ -18,7 +18,6 @@ MODULE cp_main_variables
   USE metagga,           ONLY : kedtaur, kedtaus, kedtaug
   USE atoms_type_module, ONLY : atoms_type
   USE cell_base,         ONLY : boxdimensions
-  USE charge_types,      ONLY : charge_descriptor, charge_descriptor_init
   USE wave_types,        ONLY : wave_descriptor, wave_descriptor_init
   USE energies,          ONLY : dft_energy_type
   !
@@ -80,19 +79,17 @@ MODULE cp_main_variables
   !
   ! charge densities and potentials
   !
-  REAL(DP), ALLOCATABLE :: rhoe(:,:,:,:)     ! charge density in real space
-  REAL(DP), ALLOCATABLE :: vpot(:,:,:,:)
-  TYPE (charge_descriptor)  :: desc           ! charge density descriptor
-  !
   ! rhog  = charge density in g space
   ! rhor  = charge density in r space (dense grid)
   ! rhos  = charge density in r space (smooth grid)
   ! rhopr   since rhor is overwritten in vofrho,
   !         this array is used to save rhor for restart file
+  ! vpot  = potential in r space (dense grid)
   !
   COMPLEX(DP), ALLOCATABLE :: rhog(:,:)
   REAL(DP),    ALLOCATABLE :: rhor(:,:), rhos(:,:)
   REAL(DP),    ALLOCATABLE :: rhopr(:,:)  
+  REAL(DP),    ALLOCATABLE :: vpot(:,:)
   !
   TYPE (wave_descriptor) :: wfill, wempt    ! wave function descriptor
                                             ! for filled and empty states
@@ -162,10 +159,11 @@ MODULE cp_main_variables
       !
       ALLOCATE( ema0bg( ngw ) )
       !
+      ALLOCATE( rhor( nnr, nspin ) )
+      !
       IF( program_name == 'CP90' ) THEN
          !
          ALLOCATE( rhopr( nnr,   nspin ) )
-         ALLOCATE( rhor( nnr,   nspin ) )
          ALLOCATE( rhos( nnrsx, nspin ) )
          ALLOCATE( rhog( ng,    nspin ) )
          !
@@ -179,12 +177,7 @@ MODULE cp_main_variables
          !
       ELSE IF( program_name == 'FPMD' ) THEN
          !
-         ALLOCATE( rhoe( nr1x, nr2x, npl, nspin ) )
+         ALLOCATE( vpot( nnr, nspin ) )
-         !
-         CALL charge_descriptor_init( desc, nr1, nr2, nr3, &
-             nr1, nr2, npl, nr1x, nr2x, npl, nspin )
-         !
-         ALLOCATE( vpot( nr1x, nr2x, npl, nspin ) )
          !
       END IF
       !
@@ -241,7 +234,7 @@ MODULE cp_main_variables
       IF( ALLOCATED( kedtaur ) ) DEALLOCATE( kedtaur )
       IF( ALLOCATED( kedtaus ) ) DEALLOCATE( kedtaus )
       IF( ALLOCATED( kedtaug ) ) DEALLOCATE( kedtaug )
-      IF( ALLOCATED( rhoe ) )    DEALLOCATE( rhoe )
+      ! IF( ALLOCATED( rhoe ) )    DEALLOCATE( rhoe )
       IF( ALLOCATED( vpot ) )    DEALLOCATE( vpot )
       IF( ALLOCATED( occn ) )    DEALLOCATE( occn )
       !

CPV/move_electrons.f90

@@ -41,7 +41,7 @@ SUBROUTINE move_electrons( nfi, tfirst, tlast, b1, b2, b3, fion, &
   USE runcp_module,         ONLY : runcp_uspp
   USE wave_constrains,      ONLY : interpolate_lambda
   USE gvecw,                ONLY : ngw
-  ! USE para_mod,             ONLY : 
+  USE orthogonalize_base,   ONLY : calphi
   !
   IMPLICIT NONE
   !
@@ -143,7 +143,7 @@ SUBROUTINE move_electrons( nfi, tfirst, tlast, b1, b2, b3, fion, &
      ! ... calphi calculates phi
      ! ... the electron mass rises with g**2
      !
-     CALL calphi( c0, ngw, ema0bg, bec, nkb, vkb, phi, nbsp )
+     CALL calphi( c0, ngw, bec, nkb, vkb, phi, nbsp, ema0bg )
      !
      ! ... begin try and error loop (only one step!)
      !

CPV/nl_base.f90

@@ -337,7 +337,7 @@
       use cvan,           only : ish
       use uspp_param,     only : nhm, nh
       use uspp,           only : nkb, dvan
-      use electrons_base, only : n => nbsp, nspin, ispin => fspin, f
+      use electrons_base, only : n => nbsp, nspin, ispin, f
       use ions_base,      only : nsp, nat, na
       !
       implicit none
@@ -400,7 +400,7 @@
       use uspp_param,     only : nhm, nh
       use cvan,           only : ish, nvb
       use ions_base,      only : nax, nat, nsp, na
-      use electrons_base, only : n => nbsp, ispin => fspin, f
+      use electrons_base, only : n => nbsp, ispin, f
       use gvecw,          only : ngw
 !
       implicit none
@@ -641,7 +641,7 @@ subroutine dennl( bec, denl )
   use cdvan,      ONLY : drhovan, dbec
   use ions_base,  only : nsp, na
   !
-  use electrons_base,     only : n => nbsp, ispin => fspin, f, nspin
+  use electrons_base,     only : n => nbsp, ispin, f, nspin
   use reciprocal_vectors, only : gstart
 
   implicit none
@@ -710,7 +710,7 @@ subroutine nlfq( c, eigr, bec, becdr, fion )
   use uspp_param,     only : nhm, nh
   use cvan,           only : ish, nvb
   use ions_base,      only : nax, nat, nsp, na
-  use electrons_base, only : n => nbsp, ispin => fspin, f
+  use electrons_base, only : n => nbsp, ispin, f
   use gvecw,          only : ngw
   use constants,      only : pi, fpi
   !

CPV/nlcc.f90

@@ -150,7 +150,7 @@
 
      integer :: nsp
      COMPLEX(DP) :: rhoetg(:)
-     REAL(DP)    :: rhoetr(:,:,:)
+     REAL(DP)    :: rhoetr(:)
      REAL(DP)    :: rhoc(:,:)
      COMPLEX(DP), INTENT(IN) :: sfac(:,:)
           

CPV/optical.f90

@@ -90,7 +90,7 @@
           COMPLEX(DP), INTENT(INOUT) :: ce(:,:,:,:)
           TYPE(wave_descriptor), INTENT(IN) :: wempt, wfill
           REAL(DP), INTENT(IN) :: occ(:,:,:)
-          REAL (DP), INTENT(in) ::   vpot(:,:,:,:)
+          REAL (DP), INTENT(in) ::   vpot(:,:)
           REAL (DP) ::   bec(:,:)
           COMPLEX(DP) :: eigr(:,:)
 
@@ -139,7 +139,7 @@
             CALL nlsm1 ( nb_l, 1, nspnl, eigr, cf(1,1,1,ispin), bec )
 
             CALL dforce_all( ispin, cf(:,:,1,ispin), wfill, occ(:,1,ispin), eforce(:,:,1), &
-                             vpot(:,:,:,ispin), eigr, bec )
+                             vpot(:,ispin), eigr, bec )
 
             CALL kohn_sham( ispin, cf(:,:,:,ispin), wfill, eforce )
 
@@ -156,7 +156,7 @@
             CALL nlsm1 ( nb_l, 1, nspnl, eigr, ce(1,1,1,ispin), bece )
             !
             CALL dforce_all( ispin, ce(:,:,1,ispin), wempt, ff( :, 1), eforce(:,:,1), &
-                             vpot(:,:,:,ispin), eigr, bece )
+                             vpot(:,ispin), eigr, bece )
             !
             CALL kohn_sham( ispin, ce(:,:,:,ispin), wempt, eforce )
 

CPV/potentials.f90

@@ -176,7 +176,7 @@
 !  -------------------------------------------------------------------------
 
       SUBROUTINE kspotential &
-        ( nfi, tprint, tforce, tstress, rhoe, desc, atoms, bec, becdr, eigr, &
+        ( nfi, tprint, tforce, tstress, rhoe, atoms, bec, becdr, eigr, &
           ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht, fi, vpot, edft, timepre )
 
         USE charge_density,    ONLY: rhoofr
@@ -185,7 +185,6 @@
         USE cell_module,       ONLY: boxdimensions
         USE atoms_type_module, ONLY: atoms_type
         USE wave_types,        ONLY: wave_descriptor
-        USE charge_types,      ONLY: charge_descriptor
 
 ! ...   declare subroutine arguments
         LOGICAL   :: tcel
@@ -193,7 +192,7 @@
         TYPE (atoms_type),    INTENT(INOUT) :: atoms
         COMPLEX(DP),         INTENT(INOUT) :: c0(:,:,:,:)
         TYPE (wave_descriptor),  INTENT(IN) :: cdesc
-        REAL(DP) :: rhoe(:,:,:,:)
+        REAL(DP) :: rhoe(:,:)
         COMPLEX(DP) :: ei1(:,:)
         COMPLEX(DP) :: ei2(:,:)
         COMPLEX(DP) :: ei3(:,:)
@@ -203,17 +202,16 @@
         REAL(DP) :: bec(:,:)
         REAL(DP) :: becdr(:,:,:)
         TYPE (dft_energy_type) :: edft
-        REAL(DP)    :: vpot(:,:,:,:)
+        REAL(DP)    :: vpot(:,:)
         COMPLEX(DP), INTENT(IN) :: sfac(:,:)
         LOGICAL, INTENT(IN) :: tforce, tstress, tprint
         REAL(DP), INTENT(OUT) :: timepre
-        TYPE (charge_descriptor),  INTENT(IN) :: desc
 
         edft%enl = nlrh_m( c0, cdesc, tforce, atoms, fi, bec, becdr, eigr )
 
-        CALL rhoofr( nfi, c0, cdesc, fi, rhoe, desc, ht )
+        CALL rhoofr( nfi, c0, cdesc, fi, rhoe, ht )
 
-        CALL vofrhos( tprint, tforce, tstress, rhoe, desc, atoms, vpot, bec, &
+        CALL vofrhos( tprint, tforce, tstress, rhoe, atoms, vpot, bec, &
                       c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre,  &
                       ht, edft )
 
@@ -223,7 +221,7 @@
 !=----------------------------------------------------------------------------=!
 
    SUBROUTINE vofrhos &
-      ( tprint, tforce, tstress, rhoe, desc, atoms, vpot, bec, c0, cdesc, fi, &
+      ( tprint, tforce, tstress, rhoe, atoms, vpot, bec, c0, cdesc, fi, &
         eigr, ei1, ei2, ei3, sfac, timepre, box, edft )
 
       !  this routine computes:
@@ -260,7 +258,7 @@
 
       ! ... include modules
 
-      USE control_flags,  ONLY: tscreen, tchi2, iprsta
+      USE control_flags,  ONLY: tscreen, tchi2, iprsta, force_pairing
       USE mp_global,      ONLY: nproc, mpime, root, group
       USE mp,             ONLY: mp_sum
       USE cell_module,    ONLY: boxdimensions
@@ -274,9 +272,6 @@
       USE charge_density, ONLY: gradrho
       USE chi2,           ONLY: rhochi, allocate_chi2, deallocate_chi2
       USE vanderwaals,    ONLY: tvdw, vdw
-      USE charge_density, ONLY: checkrho
-      USE charge_types,   ONLY: charge_descriptor
-      USE wave_functions, ONLY: dft_kinetic_energy
       USE wave_types,     ONLY: wave_descriptor
       USE io_global,      ONLY: ionode, stdout
       USE sic_module,     ONLY: self_interaction, sic_epsilon, sic_alpha !!TO ADD!!!
@@ -285,15 +280,17 @@
       USE atom,           ONLY: nlcc
       USE core,           ONLY: nlcc_any, rhocg
       USE core,           ONLY: add_core_charge, core_charge_forces
+      !
       USE reciprocal_vectors,        ONLY: gx
       USE atoms_type_module,         ONLY: atoms_type
       USE exchange_correlation,      ONLY: exch_corr_energy
+      use grid_dimensions,           only: nr1, nr2, nr3, nnrx
 
       IMPLICIT NONE
 
 ! ... declare subroutine arguments
       LOGICAL, INTENT(IN) :: tprint, tforce, tstress
-      REAL(DP)            :: vpot(:,:,:,:)
+      REAL(DP)            :: vpot(:,:)
       REAL(DP),    INTENT(IN) :: fi(:,:,:)
       REAL(DP)    :: bec(:,:)
       COMPLEX(DP) :: ei1(:,:)
@@ -303,10 +300,9 @@
       COMPLEX(DP),   INTENT(IN) :: c0(:,:,:,:)
       TYPE (atoms_type), INTENT(INOUT) :: atoms
       TYPE (wave_descriptor), INTENT(IN) :: cdesc
-      TYPE (charge_descriptor),    INTENT(IN) :: desc
       TYPE (boxdimensions),    INTENT(INOUT) :: box
       TYPE (dft_energy_type) :: edft
-      REAL(DP) :: rhoe(:,:,:,:)
+      REAL(DP) :: rhoe(:,:)
       COMPLEX(DP), INTENT(IN) :: sfac(:,:)
 
       TYPE (dft_energy_type) :: edft_self
@@ -321,16 +317,16 @@
       COMPLEX(DP), ALLOCATABLE :: vloc(:), self_vloc(:)
       COMPLEX(DP), ALLOCATABLE :: rho12(:), rhoeg(:,:), self_rhoeg(:)
       COMPLEX(DP), ALLOCATABLE :: rhoetg(:,:)
-      REAL(DP), ALLOCATABLE :: rhoetr(:,:,:,:)
+      REAL(DP), ALLOCATABLE :: rhoetr(:,:)
       REAL(DP), ALLOCATABLE :: fion_vdw(:,:)
-      REAL(DP), ALLOCATABLE :: grho(:,:,:,:,:)
+      REAL(DP), ALLOCATABLE :: grho(:,:,:)
-      REAL(DP), ALLOCATABLE :: v2xc(:,:,:,:,:)
+      REAL(DP), ALLOCATABLE :: v2xc(:,:,:)
       REAL(DP), ALLOCATABLE :: fion(:,:)
 
-      REAL(DP), ALLOCATABLE :: self_rho(:,:,:,:)
+      REAL(DP), ALLOCATABLE :: self_rho(:,:)
-      REAL(DP), ALLOCATABLE :: self_vpot(:,:,:,:)
+      REAL(DP), ALLOCATABLE :: self_vpot(:,:)
-      REAL(DP), ALLOCATABLE :: self_grho(:,:,:,:,:)
+      REAL(DP), ALLOCATABLE :: self_grho(:,:,:)
-      REAL(DP), ALLOCATABLE :: self_v2xc(:,:,:,:,:)
+      REAL(DP), ALLOCATABLE :: self_v2xc(:,:,:)
 
       REAL(DP) ::  pail(3,3)
 
@@ -347,35 +343,21 @@
 
       LOGICAL :: ttscreen, ttsic, tgc
 
-      INTEGER ig1, ig2, ig3, is, ia, ig, isc, iflag, ispin
+      INTEGER ig1, ig2, ig3, is, ia, ig, isc, iflag, iss
-      INTEGER ik, i, j, k, isa, idum, nspin
+      INTEGER ik, i, j, k, isa, idum, nspin, iswfc
-      INTEGER nr1_l, nr2_l, nr3_l, nr1_g, nr2_g, nr3_g,  nnr_l
-      INTEGER :: nr1x, nr2x, nr3x
       INTEGER :: ierr
 
       DATA iflag / 0 /
       SAVE iflag, desr
 
+      REAL(DP), EXTERNAL :: enkin
+
 !  end of declarations
 !  ----------------------------------------------
 
       IF(timing) s0 = cclock()
 
-      nspin = desc % nspin
+      nspin = SIZE( rhoe, 2 )
-      
-      nr1_l = desc % nxl
-      nr2_l = desc % nyl
-      nr3_l = desc % nzl
-
-      nr1_g = desc % nx
-      nr2_g = desc % ny
-      nr3_g = desc % nz
-
-      nnr_l    = nr1_l * nr2_l * nr3_l
-
-      nr1x = dfftp%nr1x
-      nr2x = dfftp%nr2x
-      nr3x = dfftp%npl
       
       edft%evdw = 0.0d0
       !
@@ -392,7 +374,9 @@
         CALL allocate_chi2(ngm)
       END IF
 
-      ALLOCATE( rhoetr( nr1x, nr2x, nr3x, nspin) )
+      ALLOCATE( rhoetr( nnrx, nspin ) )
+      rhoetr = 0.0d0
+
       ALLOCATE( fion( 3, atoms%nat ) )
 
       fion = atoms%for( 1:3, 1:atoms%nat )
@@ -400,44 +384,43 @@
       pail = box%pail
 
       IF(tgc) THEN
-        ALLOCATE( grho( nr1x, nr2x, nr3x, 3, nspin ) )
+        ALLOCATE( grho( nnrx, 3, nspin ) )
-        ALLOCATE( v2xc( nr1x, nr2x, nr3x, nspin, nspin) )
+        ALLOCATE( v2xc( nnrx, nspin, nspin) )
       ELSE
-        ALLOCATE( grho( 1, 1, 1, 1, 1 ) )
+        ALLOCATE( grho( 1, 1, 1 ) )
-        ALLOCATE( v2xc( 1, 1, 1, 1, 1 ) )
+        ALLOCATE( v2xc( 1, 1, 1 ) )
       END IF
+      grho = 0.0d0
+      v2xc = 0.0d0
 
       ALLOCATE( rhoeg(ngm, nspin) )
       ALLOCATE( rhoetg(ngm, nspin) )
 
-      !edft%self_sxc  = 0.d0
-      !edft%sxc       = 0.d0
-      !edft%self_ehte = 0.d0
-      !edft%eht       = 0.d0
-
-    
       IF( ttsic ) THEN
       
         IF ( tgc ) THEN
 
-          ALLOCATE(self_grho( nr1x, nr2x, nr3x, 3, nspin ), STAT = ierr)
+          ALLOCATE(self_grho( nnrx, 3, nspin ), STAT = ierr)
           IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_grho ', ierr)
     
-          ALLOCATE(self_v2xc( nr1x, nr2x, nr3x, nspin, nspin ), STAT = ierr)
+          ALLOCATE(self_v2xc( nnrx, nspin, nspin ), STAT = ierr)
           IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_v2xc ', ierr)
   
+          self_grho  = 0.D0
           self_v2xc  = 0.D0
 
         END IF !on tgc
 
-        ALLOCATE (self_vpot( nr1x, nr2x, nr3x, 2 ), STAT = ierr)
+        ALLOCATE (self_vpot( nnrx, 2 ), STAT = ierr)
         IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_vpot ', ierr)
    
         self_vpot  = 0.D0
      
-        ALLOCATE (self_rho( nr1x, nr2x, nr3x, 2), STAT = ierr)
+        ALLOCATE (self_rho( nnrx, 2), STAT = ierr)
         IF( ierr /= 0 ) CALL errore(' vofrhos ', ' allocating self_rho ', ierr)
 
+        self_rho  = 0.D0
+
       END IF !on self_interaction
 
       IF(timing) s1 = cclock()
@@ -446,7 +429,12 @@
 
       edft%ekin  = 0.0_DP
       edft%emkin = 0.0_DP
-      edft%ekin  = dft_kinetic_energy(c0, cdesc, fi, edft%emkin)
+
+      DO iss = 1, nspin
+        iswfc = iss
+        IF( force_pairing ) iswfc = 1
+        edft%ekin  = edft%ekin + enkin( c0(1,1,1,iswfc), SIZE(c0,1), fi(1,1,iss), cdesc%nbl(iss) )
+      END DO
 
       IF(tprint) THEN
         IF( ionode .AND.  ttscreen ) &
@@ -454,11 +442,6 @@
       END IF
 
 
-! ... reciprocal-space vectors are in units of alat/(2 pi) so a
-! ... multiplicative factor (2 pi/alat)**2 is required
-      edft%ekin  = edft%ekin  * tpiba2
-      edft%emkin = edft%emkin * tpiba2
-
       IF( tstress .OR. tforce .OR. iflag == 0 )  THEN
         CALL vofesr( edft%esr, desr, fion, atoms, tstress, box )
         IF( iflag == 0 ) &
@@ -469,23 +452,26 @@
 
       IF(timing) s2 = cclock()
 
-! ... FFT: rho(r) --> rho(g)  
+      DO iss = 1, nspin
-      DO ispin = 1, nspin
 
-        CALL pfwfft( rhoeg(:,ispin), rhoe(:,:,:,ispin) )
+        ! ... FFT: rho(r) --> rho(g)  
+
+        CALL pfwfft( rhoeg(:,iss), rhoe(:,iss) )
 
 
         ! ... add core contribution to the charge
 
-        CALL ZCOPY( SIZE(rhoetg,1), rhoeg(1,ispin), 1, rhoetg(1,ispin), 1 )
+        CALL ZCOPY( SIZE(rhoeg,1) , rhoeg(1,iss), 1, rhoetg(1,iss), 1 )
-        CALL DCOPY( SIZE(rhoe(:,:,:,ispin)), rhoe(1,1,1,ispin), 1, rhoetr(1,1,1,ispin), 1 )
+        CALL DCOPY( nnrx, rhoe(1,iss), 1, rhoetr(1,iss), 1 )
 
         IF( nlcc_any ) THEN
 
           ! ...     add core correction
           ! ...     rhoetg = rhoeg + cc
           ! ...     rhoetr = rhoe  + cc
-          CALL add_core_charge( rhoetg(:,ispin), rhoetr(:,:,:,ispin), sfac, rhocg, atoms%nsp )
+
+          CALL add_core_charge( rhoetg(:,iss), rhoetr(:,iss), sfac, rhocg, atoms%nsp )
+
         ELSE
 
           ! ...     no core correction
@@ -493,6 +479,7 @@
           ! ...     rhoetr = rhoe
 
           ! ...     chi2
+
           IF(tchi2) THEN
             IF(nspin.GT.1) CALL errore(' vofrho ',' spin + tchi ',nspin)
             rhochi = rhoeg(:,1)
@@ -501,7 +488,7 @@
         END IF
 
         IF(tgc) THEN
-          CALL gradrho( rhoetg(:,ispin), grho(:,:,:,:,ispin), gx )
+          CALL gradrho( rhoetg(:,iss), grho(:,:,iss), gx )
         END IF
 
       END DO
@@ -518,26 +505,26 @@
       !
       IF ( ttsic ) THEN                
 
-         self_rho(:,:,:,1) = rhoetr(:,:,:,2)
+         self_rho(:,1) = rhoetr(:,2)
-         self_rho(:,:,:,2) = rhoetr(:,:,:,2)
+         self_rho(:,2) = rhoetr(:,2)
 
          IF (tgc) THEN
-            self_grho(:,:,:,:,1) = grho(:,:,:,:,2)
+            self_grho(:,:,1) = grho(:,:,2)
-            self_grho(:,:,:,:,2) = grho(:,:,:,:,2)
+            self_grho(:,:,2) = grho(:,:,2)
          ENDIF
 
          CALL exch_corr_energy( self_rho, rhoetg, self_grho, self_vpot, &
                  self_sxcp, self_vxc, self_v2xc )
 
-         vpot (:,:,:,1) = ( 1.0d0 - sic_alpha ) * vpot(:,:,:,1)
+         vpot (:,1) = ( 1.0d0 - sic_alpha ) * vpot(:,1)
          ! 
-         vpot (:,:,:,2) = ( 1.0d0 - sic_alpha ) * vpot(:,:,:,2) + sic_alpha * ( self_vpot(:,:,:,2) + self_vpot(:,:,:,1) )
+         vpot (:,2) = ( 1.0d0 - sic_alpha ) * vpot(:,2) + sic_alpha * ( self_vpot(:,2) + self_vpot(:,1) )
 
       IF (tgc) THEN
         !
-        v2xc(:,:,:,1,1) = ( 1.0d0 - sic_alpha ) * v2xc(:,:,:,1,1)
+        v2xc(:,1,1) = ( 1.0d0 - sic_alpha ) * v2xc(:,1,1)
         !
-        v2xc(:,:,:,2,2) = ( 1.0d0 - sic_alpha ) * v2xc(:,:,:,2,2) +sic_alpha * ( self_v2xc(:,:,:,2,2) + self_v2xc(:,:,:,1,1) )
+        v2xc(:,2,2) = ( 1.0d0 - sic_alpha ) * v2xc(:,2,2) +sic_alpha * ( self_v2xc(:,2,2) + self_v2xc(:,1,1) )
         !
       END IF
 
@@ -552,14 +539,14 @@
       END IF  
 
       IF ( tstress ) THEN
-        strvxc = ( edft%sxc - vxc ) * omega / DBLE( nr1_g * nr2_g * nr3_g )
+        strvxc = ( edft%sxc - vxc ) * omega / DBLE( nr1 * nr2 * nr3 )
       END IF
 
       IF( nlcc_any ) THEN
 ! ...   xc potential (vpot) from real to G space, to compute nlcc forces
 ! ...   rhoetg = fwfft(vpot)
-        DO ispin = 1, nspin
+        DO iss = 1, nspin
-          CALL pfwfft( rhoetg(:,ispin), vpot(:,:,:,ispin) )
+          CALL pfwfft( rhoetg(:,iss), vpot(:,iss) )
         END DO
 ! ...   now rhoetg contains the xc potential
         IF (tforce) THEN
@@ -584,7 +571,7 @@
       !
       CALL vofloc(ttscreen, tforce, edft%ehte, edft%ehti, ehp, & 
            eps, vloc, rhoeg, fion, atoms, rhops, vps, eigr, &
-           ei1, ei2, ei3, sfac, box, desc )
+           ei1, ei2, ei3, sfac, box )
       !
       edft%self_ehte = 0.d0
 
@@ -600,16 +587,16 @@
 
         !  working on the total charge density
 
-        CALL self_vofloc( ttscreen, self_ehtep, self_vloc, self_rhoeg, box, desc)
+        CALL self_vofloc( ttscreen, self_ehtep, self_vloc, self_rhoeg, box )
         !
-        CALL pinvfft( self_vpot(:,:,:,1), self_vloc(:))
+        CALL pinvfft( self_vpot(:,1), self_vloc(:))
 
-        self_vpot(:,:,:,1) = sic_epsilon * self_vpot(:,:,:,1)
+        self_vpot(:,1) = sic_epsilon * self_vpot(:,1)
         !
         edft%self_ehte     = sic_epsilon * DBLE( self_ehtep )
  
-        vpot(:,:,:,1) =  vpot(:,:,:,1) - self_vpot(:,:,:,1)
+        vpot(:,1) =  vpot(:,1) - self_vpot(:,1)
-        vpot(:,:,:,2) =  vpot(:,:,:,2) + self_vpot(:,:,:,1)
+        vpot(:,2) =  vpot(:,2) + self_vpot(:,1)
 
         DEALLOCATE( self_vloc, self_rhoeg )
 
@@ -628,13 +615,13 @@
       IF(timing) s5 = cclock()
 
 
-      DO ispin = 1, nspin
+      DO iss = 1, nspin
 
 ! ...   add hartree end local pseudo potentials ( invfft(vloc) )
 ! ...   to xc potential (vpot).
 ! ...   vpot = vpot + invfft(vloc)
 
-        CALL pinvfft( vpot(:,:,:,ispin), vloc(:), 1.0d0 )
+        CALL pinvfft( vpot(:,iss), vloc(:), 1.0d0 )
 
       END DO
 
@@ -660,10 +647,10 @@
       CALL mp_sum(edft%ehte, group)
       CALL mp_sum(edft%ehti, group)
       CALL mp_sum(edft%self_ehte, group)
-      CALL mp_sum(edft%ekin, group)
+      ! CALL mp_sum(edft%ekin, group)  ! already summed up
       CALL mp_sum(edft%emkin, group)
 
-      CALL total_energy(edft,omega,vxc,eps,self_vxc,nr1_g*nr2_g*nr3_g)
+      CALL total_energy(edft,omega,vxc,eps,self_vxc,nr1*nr2*nr3)
 
 !fran: the output is introduced only in the print_energies.f90
 !fran: in this way you print only each print_step
@@ -761,21 +748,20 @@
 
 !=----------------------------------------------------------------------------=!
 
-  SUBROUTINE cluster_bc( screen_coul, hg, box, desc )
+  SUBROUTINE cluster_bc( screen_coul, hg, box )
 
       USE green_functions, ONLY: greenf
       USE mp_global, ONLY: mpime
       USE fft, ONLY : pfwfft
       USE fft_base, ONLY: dfftp
-      USE charge_types, ONLY: charge_descriptor
       USE processors_grid_module, ONLY: get_grid_info
       USE cell_module, ONLY: boxdimensions, s_to_r, alat
       USE constants, ONLY: gsmall, pi
       USE cell_base, ONLY: tpiba2
+      use grid_dimensions, only: nr1, nr2, nr3, nr1l, nr2l, nr3l, nnrx
       
       REAL(DP), INTENT(IN) :: hg(:)
       TYPE (boxdimensions),    INTENT(IN) :: box
-      TYPE (charge_descriptor),    INTENT(IN) :: desc
       COMPLEX(DP) :: screen_coul(:)
 
 ! ... declare external function
@@ -783,21 +769,12 @@
       EXTERNAL erf, erfc
 
 ! ... Locals
-      REAL(DP), ALLOCATABLE :: grr(:,:,:)
+      REAL(DP), ALLOCATABLE :: grr(:)
       COMPLEX(DP), ALLOCATABLE :: grg(:)
       REAL(DP) :: rc, r(3), s(3), rmod, g2, rc2, arg, omega, fact
-      INTEGER   :: ig, i, j, k
+      INTEGER   :: ig, i, j, k, ir
-      INTEGER   :: nr1_l, nr2_l, nr3_l, nr1_g, nr2_g, nr3_g
       INTEGER   :: ir1, ir2, ir3
 
-
-      nr1_l = desc % nxl
-      nr2_l = desc % nyl
-      nr3_l = desc % nzl
-      nr1_g = desc % nx
-      nr2_g = desc % ny
-      nr3_g = desc % nz
-
       ir1 = 1
       ir2 = 1
       ir3 = 1
@@ -805,28 +782,31 @@
         ir3 = ir3 + dfftp%npp( k )
       END DO
 
-      ALLOCATE( grr( dfftp%nr1x, dfftp%nr2x, dfftp%npl ) )
+      ALLOCATE( grr( nnrx ) )
       ALLOCATE( grg( SIZE( screen_coul ) ) )
 
+      grr = 0.0d0
+
 ! ... Martina and Tuckerman convergence criterium
       rc  = 7.0d0 / alat
       rc2 = rc**2
       omega = box%deth
-      fact  = omega / ( nr1_g * nr2_g * nr3_g )
+      fact  = omega / ( nr1 * nr2 * nr3 )
-      IF( MOD(nr1_g * nr2_g * nr3_g, 2) /= 0 ) fact = -fact
+      IF( MOD(nr1 * nr2 * nr3, 2) /= 0 ) fact = -fact
 
-      DO k = 1, nr3_l
+      DO k = 1, nr3l
-        s(3) = DBLE ( (k-1) + (ir3 - 1) ) / nr3_g - 0.5d0
+        s(3) = DBLE ( (k-1) + (ir3 - 1) ) / nr3 - 0.5d0
-        DO j = 1, nr2_l
+        DO j = 1, nr2l
-          s(2) = DBLE ( (j-1) + (ir2 - 1) ) / nr2_g - 0.5d0
+          s(2) = DBLE ( (j-1) + (ir2 - 1) ) / nr2 - 0.5d0
-          DO i = 1, nr1_l
+          DO i = 1, nr1l
-            s(1) = DBLE ( (i-1) + (ir1 - 1) ) / nr1_g - 0.5d0
+            s(1) = DBLE ( (i-1) + (ir1 - 1) ) / nr1 - 0.5d0
             CALL S_TO_R(S, R, box)
             rmod = SQRT( r(1)**2 + r(2)**2 + r(3)**2 )
+            ir =  i + (j-1)*dfftp%nr1x + (k-1)*dfftp%nr1x*dfftp%nr2x
             IF( rmod < gsmall ) THEN
-              grr(i,j,k) = fact * 2.0d0 * rc / SQRT( pi )
+              grr( ir ) = fact * 2.0d0 * rc / SQRT( pi )
             ELSE
-              grr(i,j,k) = fact * erf( rc * rmod ) / rmod
+              grr( ir ) = fact * erf( rc * rmod ) / rmod
             END IF
           END DO
         END DO
@@ -856,11 +836,11 @@
 !  BEGIN manual
 
       SUBROUTINE vofloc(tscreen, tforce, ehte, ehti, eh, eps, vloc, rhoeg, &
-                 fion, atoms, rhops, vps, eigr, ei1, ei2, ei3, sfac, ht, desc )
+                 fion, atoms, rhops, vps, eigr, ei1, ei2, ei3, sfac, ht )
 
 !  this routine computes:
 !  omega = ht%deth
-!  rho_e(ig)    =  (sum over ispin) rhoeg(ig,ispin) 
+!  rho_e(ig)    =  (sum over iss) rhoeg(ig,iss) 
 !  rho_I(ig)    =  (sum over is) sfac(is,ig) * rhops(ig,is) 
 !  vloc_h(ig)   =  fpi / ( g(ig) * tpiba2 ) * { rho_e(ig) + rho_I(ig) }
 !  vloc_ps(ig)  =  (sum over is) sfac(is,ig) * vps(ig,is)
@@ -884,7 +864,7 @@
 !  tx_ps(ig,is)    = vps(ig,is) * CONJG( rho_e(ig) )
 !  gx(ig)          = CMPLX(0.D0, gx(1,ig)) * tpiba
 !  fion(x,isa)     = fion(x,isa) + 
-!      Fact * omega * ( sum over ig, ispin) (tx_h(ig,is) + tx_ps(ig,is)) * 
+!      Fact * omega * ( sum over ig, iss) (tx_h(ig,is) + tx_ps(ig,is)) * 
 !      gx(ig) * eigrx(ig,isa) * eigry(ig,isa) * eigrz(ig,isa) 
 !  if Gamma symmetry Fact = 2.0 else Fact = 1
 !
@@ -895,7 +875,6 @@
       USE control_flags, ONLY: gamma_only
       USE cell_base, ONLY: tpiba2, tpiba
       USE cell_module, ONLY: boxdimensions
-      USE charge_types, ONLY: charge_descriptor
       USE atoms_type_module, ONLY: atoms_type
       USE io_global, ONLY: stdout
       USE grid_dimensions, ONLY: nr1, nr2, nr3
@@ -910,7 +889,6 @@
 
       TYPE (atoms_type) :: atoms
       TYPE (boxdimensions), INTENT(in) :: ht
-      TYPE (charge_descriptor), INTENT(IN) :: desc
       LOGICAL      :: tforce
       LOGICAL      :: tscreen
       REAL(DP)    :: fion(:,:)
@@ -927,7 +905,7 @@
 
 ! ... Locals
 
-      INTEGER      :: is, ia, isa, ig, ig1, ig2, ig3, nspin, ispin
+      INTEGER      :: is, ia, isa, ig, ig1, ig2, ig3, nspin, iss
       REAL(DP)    :: fpibg, cost, omega
       COMPLEX(DP) :: cxc, rhet, rhog, vp, rp, gxc, gyc, gzc
       COMPLEX(DP) :: teigr, cnvg, cvn, tx, ty, tz, vscreen
@@ -944,7 +922,7 @@
 
       IF( tscreen ) THEN
         ALLOCATE( screen_coul( ngm ) )
-        CALL cluster_bc( screen_coul, g, ht, desc )
+        CALL cluster_bc( screen_coul, g, ht )
       END IF
 
 !=======================================================================
@@ -1020,9 +998,9 @@
       IF(TFORCE) THEN
 ! ...   each processor add its own contribution to the array FION
         IF( gamma_only ) THEN
-          cost = 2.D0 * ht%deth * tpiba
+          cost = 2.D0 * omega * tpiba
         ELSE
-          cost = ht%deth * tpiba
+          cost = omega * tpiba
         END IF
         FION = FION + DBLE(ftmp) * cost
       END IF
@@ -1046,11 +1024,11 @@
         eh      = eh   +  vscreen *        rhog * CONJG(rhog)
         ehte    = ehte +  vscreen *   DBLE(rhet * CONJG(rhet))
         ehti    = ehti +  vscreen *   DBLE(  rp * CONJG(rp))
-        DO ispin = 1, nspin
+        DO iss = 1, nspin
           IF( gamma_only ) THEN
-            eps = eps + vp * CONJG(RHOEG(1,ispin)) * 0.5d0
+            eps = eps + vp * CONJG(RHOEG(1,iss)) * 0.5d0
           ELSE
-            eps = eps + vp * CONJG(RHOEG(1,ispin))
+            eps = eps + vp * CONJG(RHOEG(1,iss))
           END IF
         END DO
       END IF
@@ -1079,8 +1057,7 @@
       USE cell_module, ONLY: s_to_r, boxdimensions, pbcs
       USE mp_global, ONLY: nproc, mpime, group
       USE mp, ONLY: mp_sum
-      USE parallel_types, ONLY: BLOCK_PARTITION_SHAPE
+      USE parallel_types, ONLY: BLOCK_PARTITION_DIST
-      USE descriptors_module, ONLY: global_index, local_dimension
       USE atoms_type_module, ONLY: atoms_type
       USE ions_base, ONLY: rcmax, zv
  
@@ -1099,6 +1076,9 @@
       REAL(DP) :: erf, erfc
       EXTERNAL erf, erfc
 
+      INTEGER :: ldim_block, gind_block
+      EXTERNAL ldim_block, gind_block
+
       
 ! ... LOCALS 
 
@@ -1188,10 +1168,8 @@
       ESR     = 0.0_DP
       DESR    = 0.0_DP
 
-      ! NA_LOC = LOCALDIM(npt,NPROC,ME)
+      NA_LOC = ldim_block( npt, nproc, mpime)
-      NA_LOC = local_dimension( npt, 1, mpime, 0, nproc, BLOCK_PARTITION_SHAPE)
+      IA_S   = gind_block( 1, npt, nproc, mpime )
-      ! IA_S   = GLOBALINDEX(1,npt,NPROC,ME)
-      IA_S   = global_index( 1, npt, 1, mpime, 0, nproc, BLOCK_PARTITION_SHAPE )
       IA_E   = IA_S + NA_LOC - 1
 
       DO ia = ia_s, ia_e
@@ -1286,7 +1264,7 @@
 !  ----------------------------------------------
 !  BEGIN manual
 
-      SUBROUTINE self_vofloc(tscreen, ehte, vloc, rhoeg, ht, desc)
+      SUBROUTINE self_vofloc(tscreen, ehte, vloc, rhoeg, ht)
 
 !  adds the hartree part of the self interaction
 !
@@ -1296,7 +1274,6 @@
       USE constants, ONLY: fpi
       USE control_flags, ONLY: gamma_only
       USE cell_module, ONLY: boxdimensions
-      USE charge_types, ONLY: charge_descriptor
       USE cell_base, ONLY: tpiba2
       USE gvecp, ONLY: ngm
       USE reciprocal_vectors, ONLY: gstart, g
@@ -1305,7 +1282,6 @@
 
 ! ... Arguments
       TYPE (boxdimensions), INTENT(in) :: ht
-      TYPE (charge_descriptor), INTENT(IN) :: desc
       LOGICAL      :: tscreen
       COMPLEX(DP) :: vloc(:)
       COMPLEX(DP) :: rhoeg(:)
@@ -1324,7 +1300,7 @@
 
       IF( tscreen ) THEN
         ALLOCATE( screen_coul( ngm ) )
-        CALL cluster_bc( screen_coul, g, ht, desc )
+        CALL cluster_bc( screen_coul, g, ht )
       END IF
 
 !=======================================================================
@@ -1375,7 +1351,7 @@
 
 
 
-      SUBROUTINE localisation( wfc, atoms_m, ht, desc)
+      SUBROUTINE localisation( wfc, atoms_m, ht)
 
 !  adds the hartree part of the self interaction
 !
@@ -1385,7 +1361,6 @@
       USE constants, ONLY: fpi
       USE control_flags, ONLY: gamma_only
       USE cell_module, ONLY: boxdimensions, s_to_r
-      USE charge_types, ONLY: charge_descriptor
       USE atoms_type_module, ONLY: atoms_type
       USE fft, ONLY : pw_invfft, pfwfft, pinvfft
       USE sic_module, ONLY: ind_localisation, nat_localisation, print_localisation
@@ -1395,6 +1370,7 @@
       USE cell_base, ONLY: tpiba2
       USE reciprocal_vectors, ONLY: gstart, g
       USE gvecp, ONLY: ngm
+      use grid_dimensions, only: nr1, nr2, nr3, nr1l, nr2l, nr3l, nnrx
 
       IMPLICIT NONE
 
@@ -1403,7 +1379,6 @@
       COMPLEX(DP), INTENT(IN) :: wfc(:)
       TYPE (atoms_type), INTENT(in) :: atoms_m
       TYPE (boxdimensions), INTENT(in) :: ht
-      TYPE (charge_descriptor), INTENT(IN) :: desc
 
 
 ! ... Locals
@@ -1411,42 +1386,36 @@
       REAL(DP)    :: ehte
       INTEGER      :: ig, at, ia, is, isa_input, isa_sorted, isa_loc
       REAL(DP)    :: fpibg, omega, aRe, aR2, R(3)
-      INTEGER      :: Xmin, Ymin, Zmin, Xmax, Ymax, Zmax
+      INTEGER      :: Xmin, Ymin, Zmin, Xmax, Ymax, Zmax, i,j,k, ir
-      INTEGER      :: nr1_l, nr2_l, nr3_l
       REAL(DP)    :: work, work2
       COMPLEX(DP) :: rhog
-      REAL(DP), ALLOCATABLE :: density(:,:,:), psi(:,:,:)
+      REAL(DP), ALLOCATABLE :: density(:), psi(:)
-      COMPLEX(DP), ALLOCATABLE :: k_density(:), cpsi(:,:,:)
+      COMPLEX(DP), ALLOCATABLE :: k_density(:), cpsi(:)
       COMPLEX(DP) :: vscreen
       COMPLEX(DP), ALLOCATABLE :: screen_coul(:)
-      INTEGER :: nr1x, nr2x, nr3x
 
 ! ... Subroutine body ...
 
 
       IF( .FALSE. ) THEN
         ALLOCATE( screen_coul( ngm ) )
-        CALL cluster_bc( screen_coul, g, ht, desc )
+        CALL cluster_bc( screen_coul, g, ht )
       END IF
 
 
-      nr1x = dfftp%nr1x
-      nr2x = dfftp%nr2x
-      nr3x = dfftp%npl
-
       omega = ht%deth
 
-      nr1_l = desc % nxl
+      ALLOCATE( density( nnrx ) )
-      nr2_l = desc % nyl
+      ALLOCATE( psi( nnrx ) )
-      nr3_l = desc % nzl
-
-      ALLOCATE( density( nr1x, nr2x, nr3x ) )
-      ALLOCATE( psi( nr1x, nr2x, nr3x ) )
-      ALLOCATE( cpsi( nr1x, nr2x, nr3x ) )
       ALLOCATE( k_density( ngm ) )
 
-      CALL pw_invfft( cpsi(:,:,:), wfc(:), wfc(:) )
+      ALLOCATE( cpsi( nnrx ) )
+      cpsi = 0.0d0
+
+      CALL pw_invfft( cpsi(:), wfc(:), wfc(:) )
+
       psi = DBLE( cpsi )
+
       DEALLOCATE( cpsi )
 
       isa_sorted = 0
@@ -1473,34 +1442,44 @@
             !WRITE(6,*) 'ATOM ', ind_localisation( isa_input )
             !WRITE(6,*) 'POS  ', atoms_m%taus( :, isa_sorted )
 
-            work  = nr1_l
+            work  = nr1l
             work2 = sic_rloc * work
             work  = work * R(1) - work2
             Xmin  = FLOOR(work)
             work  = work + 2*work2
             Xmax  = FLOOR(work)
-            IF ( Xmax > nr1_l ) Xmax = nr1_l
+            IF ( Xmax > nr1l ) Xmax = nr1l
             IF ( Xmin < 1 ) Xmin = 1
-            work  = nr2_l
+
+            work  = nr2l
             work2 = sic_rloc * work
             work  = work * R(2) - work2
             Ymin  = FLOOR(work)
             work  = work + 2*work2
             Ymax  = FLOOR(work)
-            IF ( Ymax > nr2_l ) Ymax = nr2_l
+            IF ( Ymax > nr2l ) Ymax = nr2l
             IF ( Ymin < 1 ) Ymin = 1
-            work  = nr3_l
+
+            work  = nr3l
             work2 = sic_rloc * work
             work  = work * R(3) - work2
             Zmin  = FLOOR(work)
             work  = work + 2*work2
             Zmax  = FLOOR(work)
-            IF ( Zmax > nr3_l ) Zmax = nr3_l
+            IF ( Zmax > nr3l ) Zmax = nr3l
             IF ( Zmin < 1 ) Zmin = 1
 
             density = 0.D0
-            density( Xmin:Xmax, Ymin:Ymax, Zmin:Zmax ) = &
+
-            psi( Xmin:Xmax, Ymin:Ymax, Zmin:Zmax ) * psi( Xmin:Xmax, Ymin:Ymax, Zmin:Zmax )
+            DO k = Zmin, Zmax
+               DO j = Ymin, Ymax
+                  DO i = Xmin, Xmax
+                     ir =  i + (j-1)*dfftp%nr1x + (k-1)*dfftp%nr1x*dfftp%nr2x
+                     density( ir ) = psi( ir ) * psi( ir )
+                  END DO
+               END DO
+            END DO
+
             CALL pfwfft( k_density, density )
 
 ! ... G /= 0 elements

CPV/print_out.f90

@@ -594,17 +594,15 @@
 !=----------------------------------------------------------------------------=!
 
 
-    SUBROUTINE print_sfac( rhoe, desc, sfac )
+    SUBROUTINE print_sfac( rhoe, sfac )
 
       USE mp_global, ONLY: mpime, nproc, group
       USE mp, ONLY: mp_max, mp_get, mp_put
       USE fft, ONLY : pfwfft, pinvfft
-      USE charge_types, ONLY: charge_descriptor
       USE reciprocal_vectors, ONLY: ig_l2g, gx, g
       USE gvecp, ONLY: ngm
 
-      TYPE (charge_descriptor), INTENT(IN) :: desc
+      REAL(DP), INTENT(IN) :: rhoe(:,:)
-      REAL(DP), INTENT(IN) :: rhoe(:,:,:,:)
       COMPLEX(DP), INTENT(IN) ::  sfac(:,:)
 
       INTEGER :: nspin, ispin, ip, nsp, ngx_l, ng, is, ig
@@ -615,7 +613,7 @@
       INTEGER     , ALLOCATABLE :: ig_rcv(:)
       COMPLEX(DP), ALLOCATABLE :: sfac_rcv(:,:)
 
-        nspin = SIZE(rhoe,4)
+        nspin = SIZE(rhoe,2)
         nsp   = SIZE(sfac,2)
         ngx_l = ngm
         CALL mp_max(ngx_l, group)
@@ -627,7 +625,7 @@
         ALLOCATE(sfac_rcv(ngx_l,nsp))
 ! ...   FFT: rho(r) --> rho(g)
         DO ispin = 1, nspin
-          CALL pfwfft(rhoeg(:,ispin),rhoe(:,:,:,ispin))
+          CALL pfwfft(rhoeg(:,ispin),rhoe(:,ispin))
         END DO
         IF( ionode ) THEN
           OPEN(sfacunit, FILE=TRIM(sfac_file), STATUS='UNKNOWN')
@@ -673,20 +671,15 @@
 !=----------------------------------------------------------------------------=!
 
 
-    SUBROUTINE printacc( nfi, rhoe, desc, atoms, ht, nstep_run, avgs, avgs_run )
+    SUBROUTINE printacc( nfi, nstep_run, avgs, avgs_run )
 
       USE cell_module, ONLY: boxdimensions
       USE atoms_type_module, ONLY: atoms_type
-      USE charge_types, ONLY: charge_descriptor
 
       IMPLICIT NONE
 
       INTEGER, INTENT(IN) :: nfi, nstep_run
-      REAL(DP), intent(in) :: rhoe(:,:,:,:)
-      TYPE (charge_descriptor), intent(in) :: desc
       REAL (DP) :: avgs(:), avgs_run(:)
-      TYPE (atoms_type) :: atoms
-      TYPE (boxdimensions), intent(in) :: ht
  
       IF ( nfi < 1 ) THEN
         RETURN
@@ -744,7 +737,7 @@
    
      use kinds, only: DP
      use ensemble_dft, only: tens, ismear, z0, c0diag, becdiag, dval, zaux, e0, zx
-     use electrons_base, only: nx => nbspx, n => nbsp, ispin => fspin, f, nspin
+     use electrons_base, only: nx => nbspx, n => nbsp, ispin, f, nspin
      use electrons_base, only: nel, iupdwn, nupdwn, nudx, nelt
      use energies, only: enl, ekin
      use ions_base, only: nsp

CPV/qmatrixd.f90

@@ -33,7 +33,7 @@ subroutine qmatrixd(c0, bec0,ctable, gqq, qmat, detq)
   use reciprocal_vectors, only: ng0 => gstart
   use uspp_param, only: nh, nhm
   use uspp, only : nhsa=> nkb
-  use electrons_base, only: nx => nbspx, n => nbsp, fspin
+  use electrons_base, only: nx => nbspx, n => nbsp, ispin
   use mp, only: mp_sum
   
 
@@ -80,7 +80,7 @@ subroutine qmatrixd(c0, bec0,ctable, gqq, qmat, detq)
 ! first the local part
 
         sca=(0.,0.)
-        if(fspin(ix) == fspin(jx) ) then
+        if(ispin(ix) == ispin(jx) ) then
        
 !#ifdef NEC
 !        *vdir nodep
@@ -133,7 +133,7 @@ subroutine qmatrixd(c0, bec0,ctable, gqq, qmat, detq)
 !  now the non local vanderbilt part
             
         sca =(0.,0.)
-        if(fspin(ix)==fspin(jx)) then
+        if(ispin(ix)==ispin(jx)) then
           do is=1,nvb!loop on vanderbilt species
              do ia=1,na(is)!loop on atoms
                 do iv=1,nh(is)!loop on projectors

CPV/restart.f90

@@ -230,7 +230,7 @@
 
    SUBROUTINE writefile_fpmd( nfi, trutime, c0, cm, cdesc, occ, &
      atoms_0, atoms_m, acc, taui, cdmi, &
-     ht_m, ht_0, rho, desc, vpot)
+     ht_m, ht_0, rho, vpot)
                                                                         
         USE cell_module, only: boxdimensions, r_to_s
         USE brillouin, only: kpoints, kp
@@ -240,7 +240,6 @@
         USE atoms_type_module, ONLY: atoms_type
         USE io_global, ONLY: ionode, ionode_id
         USE io_global, ONLY: stdout
-        USE charge_types, ONLY: charge_descriptor
         USE electrons_nose, ONLY: xnhe0, xnhem, vnhe
         USE electrons_base, ONLY: nbsp, nspin
         USE cell_nose, ONLY: xnhh0, xnhhm, vnhh
@@ -257,17 +256,15 @@
         REAL(DP), INTENT(IN) :: occ(:,:,:)
         TYPE (boxdimensions), INTENT(IN) :: ht_m, ht_0
         TYPE (atoms_type), INTENT(IN) :: atoms_0, atoms_m
-        REAL(DP), INTENT(IN) :: rho(:,:,:,:)
+        REAL(DP), INTENT(IN) :: rho(:,:)
-        TYPE (charge_descriptor), INTENT(IN) :: desc
         TYPE (wave_descriptor) :: cdesc
-        REAL(DP), INTENT(INOUT) :: vpot(:,:,:,:)
+        REAL(DP), INTENT(INOUT) :: vpot(:,:)
                                                                         
         REAL(DP), INTENT(IN) :: taui(:,:)
         REAL(DP), INTENT(IN) :: acc(:), cdmi(:) 
         REAL(DP), INTENT(IN) :: trutime
 
         REAL(DP), ALLOCATABLE :: lambda(:,:)
-        REAL(DP), ALLOCATABLE :: rhow(:,:)
         REAL(DP) :: ekincm
         INTEGER   :: i, j, k, iss, ir
              
@@ -278,36 +275,16 @@
         !   properties on the writefile subroutine
 
         ALLOCATE( lambda(nbsp,nbsp) )
-        ALLOCATE( rhow( nr1x * nr2x * SIZE( rho, 3 ), nspin ) )
         lambda  = 0.0d0
         ekincm = 0.0d0
         !
         !
-        IF( SIZE( rho, 1 ) /= nr1x .OR.  SIZE( rho, 2 ) /= nr2x ) THEN
-           WRITE( stdout, * ) nr1x, nr2x
-           WRITE( stdout, * ) SIZE( rho, 1 ), SIZE( rho, 2 )
-           CALL errore( ' writefile_fpmd ', ' rho dimensions do not correspond ', 1 )
-        END IF
-        !
-        DO iss = 1, nspin
-           ir = 0
-           DO k = 1, SIZE( rho, 3 )
-              DO j = 1, nr2x
-                 DO i = 1, nr1x
-                   ir = ir + 1
-                   rhow( ir, iss ) = rho( i, j, k, iss )
-                 END DO
-              END DO
-           END DO
-        END DO
-
         CALL cp_writefile( ndw, scradir, .TRUE., nfi, trutime, acc, kp%nkpt, kp%xk, kp%weight, &
           ht_0%a, ht_m%a, ht_0%hvel, ht_0%gvel, xnhh0, xnhhm, vnhh, taui, cdmi, &
           atoms_0%taus, atoms_0%vels, atoms_m%taus, atoms_m%vels, atoms_0%for, vnhp, &
           xnhp0, xnhpm, nhpcl, nhpdim, occ, occ, lambda, lambda,  &
-          xnhe0, xnhem, vnhe, ekincm, ei, rhow, c04 = c0, cm4 = cm )
+          xnhe0, xnhem, vnhe, ekincm, ei, rho, c04 = c0, cm4 = cm )
 
-        DEALLOCATE( rhow )
         DEALLOCATE( lambda )
 
      RETURN 
@@ -319,7 +296,7 @@
 
         SUBROUTINE readfile_fpmd( nfi, trutime, &
           c0, cm, cdesc, occ, atoms_0, atoms_m, acc, taui, cdmi, &
-          ht_m, ht_0, rho, desc, vpot )
+          ht_m, ht_0, rho, vpot )
                                                                         
         use electrons_base, only: nbsp
         USE cell_module, only: boxdimensions, cell_init, r_to_s, s_to_r
@@ -336,7 +313,6 @@
         USE gvecw, ONLY: ecutwfc => ecutw
         USE gvecp, ONLY: ecutrho => ecutp
         USE fft, ONLY : pfwfft, pinvfft
-        USE charge_types, ONLY: charge_descriptor
         USE ions_base, ONLY: nat, nsp, na
         USE electrons_module, ONLY: nspin
         USE control_flags, ONLY: twfcollect, force_pairing
@@ -354,10 +330,9 @@
         REAL(DP), INTENT(INOUT) :: occ(:,:,:)
         TYPE (boxdimensions), INTENT(INOUT) :: ht_m, ht_0
         TYPE (atoms_type), INTENT(INOUT) :: atoms_0, atoms_m
-        REAL(DP), INTENT(INOUT) :: rho(:,:,:,:)
+        REAL(DP), INTENT(INOUT) :: rho(:,:)
-        TYPE (charge_descriptor), INTENT(IN) :: desc
         TYPE (wave_descriptor) :: cdesc
-        REAL(DP), INTENT(INOUT) :: vpot(:,:,:,:)
+        REAL(DP), INTENT(INOUT) :: vpot(:,:)
                                                                         
         REAL(DP), INTENT(OUT) :: taui(:,:)
         REAL(DP), INTENT(OUT) :: acc(:), cdmi(:) 

CPV/restart_sub.f90

@@ -72,6 +72,8 @@ MODULE from_restart_module
     USE cell_nose,            ONLY : xnhh0, xnhhm, vnhh, cell_nosezero
     USE phase_factors_module, ONLY : strucf
     USE cg_module,            ONLY : tcg
+    USE orthogonalize,        ONLY : ortho
+    USE orthogonalize_base,   ONLY : updatc, calphi
     !
     COMPLEX(DP) :: eigr(:,:), ei1(:,:), ei2(:,:), ei3(:,:)
     COMPLEX(DP) :: eigrb(:,:)
@@ -239,7 +241,7 @@ MODULE from_restart_module
        !
        ! ... calphi calculates phi; the electron mass rises with g**2
        !
-       CALL calphi( c0, ngw, ema0bg, bec, nkb, vkb, phi, nbsp )
+       CALL calphi( c0, ngw, bec, nkb, vkb, phi, nbsp, ema0bg )
        !
        ! ... begin try and error loop ( only one step! )
        !
@@ -252,10 +254,11 @@ MODULE from_restart_module
        !
        IF ( tortho ) THEN
           !
-          CALL ortho( eigr, cm, phi, lambda, bigr, iter, &
+          CALL ortho( eigr, cm(:,:,1,1), phi(:,:,1,1), lambda, bigr, iter, &
-                      dt2bye, ortho_eps, ortho_max, delt0, bephi, becp )
+                      dt2bye, bephi, becp )
           !
-          CALL updatc( dt2bye, lambda, phi, bephi, becp, bec, cm )
+          CALL updatc( dt2bye, nbsp, lambda, SIZE(lambda,1), phi, SIZE(phi,1), &
+                       bephi, SIZE(bephi,1), becp, bec, cm )
           !
        ELSE
           !
@@ -414,7 +417,7 @@ MODULE from_restart_module
   END SUBROUTINE from_restart_sm
   !
   !--------------------------------------------------------------------------
-  SUBROUTINE from_restart_fpmd( nfi, acc, rhoe, desc, cm, c0, cdesc, &
+  SUBROUTINE from_restart_fpmd( nfi, acc, rhoe, cm, c0, cdesc, &
                                 eigr, ei1, ei2, ei3, sfac, fi, ht_m, ht_0, &
                                 atoms_m, atoms_0, bec, becdr, vpot, edft )
     !--------------------------------------------------------------------------
@@ -454,7 +457,6 @@ MODULE from_restart_module
                                       tprnfor, tpre
     USE parameters,            ONLY : nacx
     USE atoms_type_module,     ONLY : atoms_type
-    USE charge_types,          ONLY : charge_descriptor
     USE ions_base,             ONLY : vel_srt, tau_units
     USE runcp_module,          ONLY : runcp_ncpp
     USE grid_dimensions,       ONLY : nr1, nr2, nr3
@@ -476,12 +478,11 @@ MODULE from_restart_module
     COMPLEX(DP), INTENT(INOUT) :: cm(:,:,:,:), c0(:,:,:,:)
     REAL(DP)                   :: fi(:,:,:)
     TYPE(boxdimensions)        :: ht_m, ht_0
-    REAL(DP)                   :: rhoe(:,:,:,:)
+    REAL(DP)                   :: rhoe(:,:)
-    TYPE(charge_descriptor)    :: desc
     TYPE(wave_descriptor)      :: cdesc
     REAL(DP)                   :: bec(:,:)
     REAL(DP)                   :: becdr(:,:,:)
-    REAL(DP)                   :: vpot(:,:,:,:)
+    REAL(DP)                   :: vpot(:,:)
     TYPE(dft_energy_type)      :: edft
     !
     INTEGER     :: ig, ib, i, j, k, ik, nb, is, ia, ierr, isa, iss
@@ -638,9 +639,9 @@ MODULE from_restart_module
        ! 
        edft%enl = nlrh_m( c0, cdesc, ttforce, atoms_0, fi, bec, becdr, eigr )
        !
-       CALL rhoofr( nfi, c0, cdesc, fi, rhoe, desc, ht_0 )
+       CALL rhoofr( nfi, c0, cdesc, fi, rhoe, ht_0 )
        !
-       CALL vofrhos( .true. , ttforce, tstress, rhoe, desc, &
+       CALL vofrhos( .true. , ttforce, tstress, rhoe, &
                      atoms_0, vpot, bec, c0, cdesc, fi, eigr, &
                      ei1, ei2, ei3, sfac, timepre, ht_0, edft )
        !

CPV/rsmesh.f90

@@ -18,9 +18,8 @@
 !=----------------------------------------------------------------------------=!
 
         USE kinds
-        USE parallel_types, ONLY: processors_grid, descriptor, BLOCK_PARTITION_SHAPE
+        USE parallel_types, ONLY: processors_grid, descriptor, BLOCK_PARTITION_DIST
         USE processors_grid_module, ONLY: grid_init, get_grid_info, calculate_grid_dims
-        USE descriptors_module, ONLY: desc_init, get_local_dims
         USE grid_dimensions, ONLY: nr1,  nr2,  nr3, nr1x, nr2x, nr3x
         USE grid_dimensions, ONLY: nr1l, nr2l, nr3l, nnrx
         USE smooth_grid_dimensions, ONLY: nr1s,  nr2s,  nr3s, nr1sx, nr2sx, nr3sx

CPV/runcg.f90

@@ -53,7 +53,7 @@
 !  -----------------------------------------------------------------------
 !  BEGIN manual
 
-   SUBROUTINE runcg_new(tortho, tprint, rhoe, desc, atoms_0, &
+   SUBROUTINE runcg_new(tortho, tprint, rhoe, atoms_0, &
                 bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht0, occ, ei, &
                 vpot, doions, edft, maxnstep, cgthr, tconv )
 
@@ -66,7 +66,6 @@
       USE energies, ONLY: dft_energy_type, print_energies
       USE electrons_module, ONLY: pmss, eigs, nb_l
       USE cp_electronic_mass, ONLY: emass
-      USE descriptors_module, ONLY: get_local_dims, owner_of, local_index
       USE wave_functions, ONLY: cp_kinetic_energy, proj, fixwave
       USE wave_base, ONLY: dotp, hpsi
       USE wave_constrains, ONLY: update_lambda
@@ -80,7 +79,6 @@
       USE potentials, ONLY: kspotential
       USE time_step, ONLY: delt
       USE atoms_type_module, ONLY: atoms_type
-      USE charge_types, ONLY: charge_descriptor
       USE control_flags, ONLY: force_pairing
       USE environment, ONLY: start_cclock_val
       USE reciprocal_space_mesh, ONLY: gkmask_l
@@ -93,8 +91,7 @@
       TYPE (atoms_type) :: atoms_0
       COMPLEX(DP), INTENT(INOUT) :: c0(:,:,:,:), cm(:,:,:,:), cp(:,:,:,:)
       TYPE (wave_descriptor) :: cdesc
-      TYPE (charge_descriptor) :: desc
+      REAL(DP) :: rhoe(:,:)
-      REAL(DP) :: rhoe(:,:,:,:)
       COMPLEX(DP) :: eigr(:,:)
       COMPLEX(DP) :: ei1(:,:)
       COMPLEX(DP) :: ei2(:,:)
@@ -109,7 +106,7 @@
       REAL(DP) :: cgthr
 
       REAL(DP)    :: ei(:,:,:)
-      REAL(DP)    :: vpot(:,:,:,:)
+      REAL(DP)    :: vpot(:,:)
 
 ! ... declare other variables
       LOGICAL :: ttsde, ttprint, ttforce, ttstress, gzero
@@ -178,7 +175,7 @@
 
         s1 = cclock()
 
-        CALL kspotential( 1, ttprint, ttforce, ttstress, rhoe, desc, &
+        CALL kspotential( 1, ttprint, ttforce, ttstress, rhoe, &
           atoms_0, bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, vpot, edft, timepre )
 
         s2 = cclock()
@@ -189,7 +186,7 @@
 ! ...     |d H / dPsi_j > = H |Psi_j> - Sum{i} <Psi_i|H|Psi_j> |Psi_i>
 
           CALL dforce_all( ispin, c0(:,:,1,ispin), cdesc, occ(:,1,ispin), cp(:,:,1,ispin), &
-            vpot(:,:,:,ispin), eigr, bec )
+            vpot(:,ispin), eigr, bec )
  
 ! ...     Project the gradient
           IF( gamma_symmetry ) THEN
@@ -237,7 +234,7 @@
 
         !  perform line minimization in the direction of "hacca"
 
-        CALL CGLINMIN(emin, demin, tbad, edft, cp, c0, cdesc, occ, vpot, rhoe, desc, hacca, &
+        CALL CGLINMIN(emin, demin, tbad, edft, cp, c0, cdesc, occ, vpot, rhoe, hacca, &
           atoms_0, ht0, bec, becdr, eigr, ei1, ei2, ei3, sfac)
 
         ! CALL print_energies( edft )
@@ -323,7 +320,7 @@
         DO ispin = 1, nspin
 
           CALL dforce_all( ispin, c0(:,:,1,ispin), cdesc, occ(:,1,ispin), hacca(:,:,1,ispin), &
-            vpot(:,:,:,ispin), eigr, bec )
+            vpot(:,ispin), eigr, bec )
 
           nb_g( ispin ) = cdesc%nbt( ispin )
 
@@ -363,7 +360,7 @@
 !
 ! ---------------------------------------------------------------------- !
 
-    SUBROUTINE cglinmin(emin, ediff, tbad, edft, cp, c, cdesc, occ, vpot, rhoe, desc, hacca, &
+    SUBROUTINE cglinmin(emin, ediff, tbad, edft, cp, c, cdesc, occ, vpot, rhoe, hacca, &
         atoms, ht, bec, becdr, eigr, ei1, ei2, ei3, sfac)
 
 ! ... declare modules
@@ -376,7 +373,6 @@
         USE cell_module, ONLY: boxdimensions
         USE potentials, ONLY: kspotential
         USE atoms_type_module, ONLY: atoms_type
-        USE charge_types, ONLY: charge_descriptor
         USE reciprocal_space_mesh, ONLY: gkmask_l
         USE uspp,             ONLY : vkb, nkb
 
@@ -389,8 +385,7 @@
         COMPLEX(DP), INTENT(IN) :: c(:,:,:,:)
         COMPLEX(DP), INTENT(INOUT) :: cp(:,:,:,:)
         TYPE (wave_descriptor), INTENT(IN) :: cdesc
-        TYPE (charge_descriptor) :: desc
+        REAL(DP) :: rhoe(:,:)
-        REAL(DP) :: rhoe(:,:,:,:)
         COMPLEX(DP) :: sfac(:,:)
         COMPLEX(DP) :: eigr(:,:)
         COMPLEX(DP) :: ei1(:,:)
@@ -402,7 +397,7 @@
         REAL(DP) :: becdr(:,:,:)
         TYPE (dft_energy_type) :: edft
         COMPLEX (DP) ::  hacca(:,:,:,:)
-        REAL (DP), INTENT(in) ::  vpot(:,:,:,:)
+        REAL (DP), INTENT(in) ::  vpot(:,:)
 
 !
 ! ... LOCALS
@@ -605,7 +600,7 @@
         END DO
 
 
-        CALL kspotential( 1, ttprint, ttforce, ttstress, rhoe, desc, &
+        CALL kspotential( 1, ttprint, ttforce, ttstress, rhoe, &
             atoms, bec, becdr, eigr, ei1, ei2, ei3, sfac, cp, cdesc, tcel, ht, occ, vpot, edft, timepre )
 
         cgenergy = edft%etot

CPV/runcg_ion.f90

@@ -29,7 +29,7 @@
 !  -----------------------------------------------------------------------
 !  BEGIN manual
 
-   SUBROUTINE runcg_ion(nfi, tortho, tprint, rhoe, desc, atomsp, atoms0, atomsm, &
+   SUBROUTINE runcg_ion(nfi, tortho, tprint, rhoe, atomsp, atoms0, atomsm, &
       bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
       vpot, doions, edft, etol, ftol, maxiter, sdthr, maxnstep )
 
@@ -52,7 +52,6 @@
       USE atoms_type_module, ONLY: atoms_type
       USE print_out_module
       USE parameters, ONLY: nacx
-      USE charge_types, ONLY: charge_descriptor
       USE runsd_module, ONLY: runsd
 
       IMPLICIT NONE
@@ -66,8 +65,7 @@
       TYPE (atoms_type) :: atomsm
       COMPLEX(DP), INTENT(INOUT) :: c0(:,:,:,:), cm(:,:,:,:), cp(:,:,:,:)
       TYPE (wave_descriptor) :: cdesc
-      TYPE (charge_descriptor) :: desc
+      REAL(DP) :: rhoe(:,:)
-      REAL(DP) :: rhoe(:,:,:,:)
       REAL(DP) :: bec(:,:)
       REAL(DP) :: becdr(:,:,:)
       COMPLEX(DP) :: eigr(:,:)
@@ -80,7 +78,7 @@
       TYPE (dft_energy_type) :: edft
 
       REAL(DP)    :: ei(:,:,:)
-      REAL(DP)    :: vpot(:,:,:,:)
+      REAL(DP)    :: vpot(:,:)
 
       INTEGER, INTENT(IN) :: maxnstep, maxiter
       REAL(DP), INTENT(IN) :: sdthr, etol, ftol
@@ -156,7 +154,7 @@
       s1 = cclock()
       old_clock_value = s1
 
-      CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atoms0, &
+      CALL runsd(ttortho, ttprint, ttforce, rhoe, atoms0, &
          bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
          vpot, doions, edft, maxnstep, sdthr )
 
@@ -184,7 +182,7 @@
         IF(ionode) &
           WRITE( stdout,fmt="(/,8X,'cgion: iter',I5,' line minimization along gradient starting')") iter
 
-        CALL CGLINMIN(fret, edft, cp, c0, cm, cdesc, occ, ei, vpot, rhoe, desc, xi, atomsp, atoms0, &
+        CALL CGLINMIN(fret, edft, cp, c0, cm, cdesc, occ, ei, vpot, rhoe, xi, atomsp, atoms0, &
           ht, bec, becdr, eigr, ei1, ei2, ei3, sfac, maxnstep, sdthr, displ)
 
         IF( tbad ) THEN
@@ -199,7 +197,7 @@
           IF( ionode ) WRITE( stdout, fmt='(8X,"cgion: bad step")')  ! perform steepest descent
           displ = displ / 2.0d0
 
-          CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atoms0, &
+          CALL runsd(ttortho, ttprint, ttforce, rhoe, atoms0, &
             bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
             vpot, doions, edft, maxnstep, sdthr )
         
@@ -289,7 +287,7 @@
 ! ---------------------------------------------------------------------- !
 
       SUBROUTINE cglinmin(emin, edft, cp, c0, cm, cdesc, occ, ei, vpot, &
-        rhoe, desc, hacca, atomsp, atoms0, ht, bec, becdr, eigr, ei1, ei2, ei3, sfac, &
+        rhoe, hacca, atomsp, atoms0, ht, bec, becdr, eigr, ei1, ei2, ei3, sfac, &
         maxnstep, sdthr, displ)
 
 ! ... declare modules
@@ -302,7 +300,6 @@
         USE cell_module, ONLY: boxdimensions, r_to_s
         USE atoms_type_module, ONLY: atoms_type
         USE check_stop, ONLY: check_stop_now
-        USE charge_types, ONLY: charge_descriptor
         USE runsd_module, ONLY: runsd
 
         IMPLICIT NONE
@@ -315,8 +312,7 @@
         COMPLEX(DP), INTENT(INOUT) :: cp(:,:,:,:)
         COMPLEX(DP), INTENT(INOUT) :: cm(:,:,:,:)
         TYPE (wave_descriptor) :: cdesc
-        TYPE (charge_descriptor) :: desc
+        REAL(DP) :: rhoe(:,:)
-        REAL(DP) :: rhoe(:,:,:,:)
         COMPLEX(DP) :: eigr(:,:)
         COMPLEX(DP) :: ei1(:,:)
         COMPLEX(DP) :: ei2(:,:)
@@ -326,7 +322,7 @@
         REAL(DP)  :: occ(:,:,:)
         TYPE (dft_energy_type) :: edft
         REAL (DP) ::  hacca(:,:)
-        REAL (DP), INTENT(in) ::  vpot(:,:,:,:)
+        REAL (DP), INTENT(in) ::  vpot(:,:)
         REAL(DP) :: bec(:,:)
         REAL(DP) :: becdr(:,:,:)
 
@@ -546,7 +542,7 @@
          ! ...  Calculate Forces (fion) and DFT Total Energy (edft) for the new ionic
          ! ...  positions (atomsp)
 
-           CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atomsp, &
+           CALL runsd(ttortho, ttprint, ttforce, rhoe, atomsp, &
              bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
              vpot, doions, edft, maxnstep, sdthr )
 

CPV/runcp.f90

@@ -44,7 +44,6 @@
       USE mp, ONLY: mp_sum
       USE electrons_module, ONLY:  pmss, eigs, nb_l
       USE cp_electronic_mass, ONLY: emass
-      USE descriptors_module, ONLY: get_local_dims, owner_of, local_index
       USE wave_functions, ONLY : cp_kinetic_energy
       USE wave_base, ONLY: hpsi
       USE cell_module, ONLY: boxdimensions
@@ -68,7 +67,7 @@
       REAL(DP), INTENT(IN)  ::  fi(:,:,:)
       REAL(DP), INTENT(IN)  ::  bec(:,:)
       TYPE (boxdimensions), INTENT(IN)  ::  ht
-      REAL (DP) ::  vpot(:,:,:,:)
+      REAL (DP) ::  vpot(:,:)
       REAL(DP) :: ei(:,:,:)
       REAL(DP) :: timerd, timeorto
       REAL(DP) :: ekinc(:)
@@ -191,7 +190,7 @@
       TYPE (wave_descriptor), INTENT(IN) :: cdesc
       COMPLEX(DP) :: eigr(:,:)
       REAL(DP), INTENT(IN)  ::  fi(:,:,:)
-      REAL (DP) ::  vpot(:,:,:,:)
+      REAL (DP) ::  vpot(:,:)
       REAL (DP), INTENT(IN) ::  bec(:,:)
       REAL(DP), INTENT(IN) :: fccc
       LOGICAL, OPTIONAL, INTENT(IN) :: lambda, fromscra, diis, restart
@@ -255,7 +254,7 @@
 
           DO i = 1, nb, 2
 
-            CALL dforce( i, is, c0(:,:,1,is), cdesc, fi(:,1,is), c2, c3, vpot(:,:,:,is), eigr, bec )
+            CALL dforce( i, is, c0(:,:,1,is), cdesc, fi(:,1,is), c2, c3, vpot(:,is), eigr, bec )
 
             IF( tlam ) THEN
                CALL update_lambda( i, gam( :, :,is), c0(:,:,1,is), cdesc, c2 )
@@ -286,7 +285,7 @@
 
             nb = nx
 
-            CALL dforce( nx, is, c0(:,:,1,is), cdesc, fi(:,1,is), c2, c3, vpot(:,:,:,is), eigr, bec )
+            CALL dforce( nx, is, c0(:,:,1,is), cdesc, fi(:,1,is), c2, c3, vpot(:,is), eigr, bec )
 
             IF( tlam ) THEN
                CALL update_lambda( nb, gam( :, :,is), c0(:,:,1,is), cdesc, c2 )
@@ -337,7 +336,6 @@
       USE mp, ONLY: mp_sum
       USE electrons_module, ONLY: pmss, eigs, nb_l, nupdwn, nspin
       USE cp_electronic_mass, ONLY: emass
-      USE descriptors_module, ONLY: get_local_dims, owner_of, local_index
       USE wave_functions, ONLY : cp_kinetic_energy
       USE wave_base, ONLY: wave_steepest, wave_verlet
       USE wave_base, ONLY: hpsi
@@ -362,7 +360,7 @@
       COMPLEX(DP)  ::  eigr(:,:)
       REAL(DP), INTENT(INOUT) ::  fi(:,:,:)
       TYPE (boxdimensions), INTENT(IN)  ::  ht
-      REAL (DP) ::  vpot(:,:,:,:)
+      REAL (DP) ::  vpot(:,:)
       REAL(DP) :: ei(:,:,:)
       REAL(DP), INTENT(IN) :: bec(:,:)
       REAL(DP) :: timerd, timeorto
@@ -457,8 +455,8 @@
 
         DO i = 1, nb, 2
           !
-          CALL dforce( i, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,:,:,1), eigr, bec )
+          CALL dforce( i, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,1), eigr, bec )
-          CALL dforce( i, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c4, c5, vpot(:,:,:,2), eigr, bec )
+          CALL dforce( i, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c4, c5, vpot(:,2), eigr, bec )
           !
           c2 = occup(i  , ik)* (c2 + c4)
           c3 = occup(i+1, ik)* (c3 + c5)
@@ -490,8 +488,8 @@
           !
           nb = n_unp - 1
           !
-          CALL dforce( nb, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,:,:,1), eigr, bec )
+          CALL dforce( nb, 2, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,1), eigr, bec )
-          CALL dforce( nb, 2, c0(:,:,1,1), cdesc, fi(:,1,2), c4, c5, vpot(:,:,:,2), eigr, bec )
+          CALL dforce( nb, 2, c0(:,:,1,1), cdesc, fi(:,1,2), c4, c5, vpot(:,2), eigr, bec )
 
           c2 = occup(nb , ik)* (c2 + c4)
 
@@ -509,7 +507,7 @@
         END IF
 
         !
-        CALL dforce( n_unp, 1, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,:,:,1), eigr, bec )
+        CALL dforce( n_unp, 1, c0(:,:,1,1), cdesc, fi(:,1,1), c2, c3, vpot(:,1), eigr, bec )
 
         intermed  = -2.d0 * sum( c2 * conjg( c0(:, n_unp, ik, 1 ) ) )
         intermed3 = sum(c0(:,n_unp, ik, 1) * conjg( c0(:, n_unp, ik, 1)))

CPV/rundiis.f90

@@ -27,7 +27,7 @@
 !  ----------------------------------------------
 !  BEGIN manual
 
-      SUBROUTINE rundiis(tprint, rhoe, desc, atoms, &
+      SUBROUTINE rundiis(tprint, rhoe, atoms, &
                  bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cgrad, cdesc, tcel, ht0, fi, eig, &
                  vpot, doions, edft )
 
@@ -94,7 +94,6 @@
       USE brillouin, ONLY: kpoints, kp
       USE wave_types, ONLY: wave_descriptor
       USE atoms_type_module, ONLY: atoms_type
-      USE charge_types, ONLY: charge_descriptor
       USE control_flags, ONLY: force_pairing
       use grid_dimensions,    only: nr1, nr2, nr3
       USE reciprocal_vectors, ONLY: mill_l
@@ -109,7 +108,7 @@
       TYPE (atoms_type) :: atoms
       COMPLEX(DP), INTENT(INOUT) :: c0(:,:,:,:), cm(:,:,:,:), cgrad(:,:,:,:)
       TYPE (wave_descriptor) :: cdesc
-      REAL(DP) :: rhoe(:,:,:,:)
+      REAL(DP) :: rhoe(:,:)
       REAL(DP) :: bec(:,:)
       REAL(DP) :: becdr(:,:,:)
       COMPLEX(DP) :: sfac(:,:)
@@ -117,13 +116,12 @@
       COMPLEX(DP) :: ei1(:,:)
       COMPLEX(DP) :: ei2(:,:)
       COMPLEX(DP) :: ei3(:,:)
-      TYPE (charge_descriptor) :: desc
       TYPE (boxdimensions), INTENT(INOUT) ::  ht0
       REAL(DP)  :: fi(:,:,:)
       TYPE (dft_energy_type) :: edft
 
       REAL(DP)    :: eig(:,:,:)
-      REAL(DP)    :: vpot(:,:,:,:) 
+      REAL(DP)    :: vpot(:,:) 
 
 ! ... declare other variables
       INTEGER ig, ib, j, k, ik, ngw, i, is, nrt, istate, nrl, ndiis, nowv
@@ -197,8 +195,8 @@
 
       CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms%taus, nr1, nr2, nr3, atoms%nat )
       CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
-      CALL rhoofr( 1, c0, cdesc, fi, rhoe, desc, ht0)
+      CALL rhoofr( 1, c0, cdesc, fi, rhoe, ht0)
-      CALL newrho(rhoe(:,:,:,1), drho, 0)  ! memorize density
+      CALL newrho(rhoe(:,1), drho, 0)  ! memorize density
       CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms%taus, nr1, nr2, nr3, atoms%nat )
       CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
       CALL guessc0( .NOT. kp%gamma_only, bec, c0, cm, cdesc)
@@ -236,12 +234,12 @@
 
 ! ...     self consistent energy
           edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, bec, becdr, eigr)
-          CALL rhoofr( 1, c0, cdesc, fi, rhoe, desc, ht0)
+          CALL rhoofr( 1, c0, cdesc, fi, rhoe, ht0)
-          CALL vofrhos(.FALSE., tforce, tstress, rhoe, desc, atoms, &
+          CALL vofrhos(.FALSE., tforce, tstress, rhoe, atoms, &
             vpot, bec, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht0, edft)
 
 ! ...     density upgrade
-          CALL newrho(rhoe(:,:,:,1), drho, idiis)
+          CALL newrho(rhoe(:,1), drho, idiis)
           IF (ionode) WRITE( stdout,45) idiis, edft%etot, drho
           dene   = abs(edft%etot - etot_m)
           etot_m = edft%etot
@@ -250,7 +248,7 @@
 
 ! ...     recalculate potential
           edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, bec, becdr, eigr)
-          CALL vofrhos(.FALSE., tforce, tstress, rhoe, desc, atoms, &
+          CALL vofrhos(.FALSE., tforce, tstress, rhoe, atoms, &
             vpot, bec, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht0, edft)
 
           IF( idiis /= 1 )THEN
@@ -269,7 +267,7 @@
 
           edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, bec, becdr, eigr)
 
-          CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,:,:,1), eigr, bec )
+          CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec )
 
           IF(.NOT.kp%gamma_only) THEN
             DO ik = 1, kp%nkpt
@@ -285,7 +283,7 @@
           call entropy_s(fi(1,1,1),temp_elec,cdesc%nbl(1),edft%ent)
 
           edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, bec, becdr, eigr)
-          CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,:,:,1), eigr, bec )
+          CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec )
 
           DO ik = 1, kp%nkpt
             DO ib = 1, cdesc%nbl( 1 )
@@ -298,7 +296,7 @@
 ! ...     DIIS on c0 at FIXED potential
           edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, bec, becdr, eigr)
 
-          CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,:,:,1), eigr, bec )
+          CALL dforce_all( 1, c0(:,:,1,1), cdesc, fi(:,1,1), cgrad(:,:,1,1), vpot(:,1), eigr, bec )
 
           IF( kp%gamma_only ) THEN
             CALL proj( 1, cgrad(:,:,1,1), cdesc, c0(:,:,1,1), cdesc, lambda)
@@ -388,7 +386,7 @@
 !  ----------------------------------------------
 !  BEGIN manual
 
-      SUBROUTINE runsdiis(tprint, rhoe, desc, atoms, &
+      SUBROUTINE runsdiis(tprint, rhoe, atoms, &
                  bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cgrad, cdesc, tcel, ht0, fi, eig, &
                  vpot, doions, edft )
 
@@ -451,7 +449,6 @@
       USE brillouin, ONLY: kpoints, kp
       USE wave_types
       USE atoms_type_module, ONLY: atoms_type
-      USE charge_types, ONLY: charge_descriptor
       USE local_pseudo, ONLY: vps
       USE uspp,             ONLY : vkb, nkb
 
@@ -462,8 +459,7 @@
       TYPE (atoms_type) :: atoms
       COMPLEX(DP), INTENT(INOUT) :: c0(:,:,:,:), cm(:,:,:,:), cgrad(:,:,:,:)
       TYPE (wave_descriptor) :: cdesc
-      REAL(DP) :: rhoe(:,:,:,:)
+      REAL(DP) :: rhoe(:,:)
-      TYPE (charge_descriptor) :: desc
       COMPLEX(DP) :: eigr(:,:)
       COMPLEX(DP) :: ei1(:,:)
       COMPLEX(DP) :: ei2(:,:)
@@ -476,7 +472,7 @@
       TYPE (dft_energy_type) :: edft
 
       REAL(DP)    :: eig(:,:,:)
-      REAL(DP)    :: vpot(:,:,:,:)
+      REAL(DP)    :: vpot(:,:)
 
 ! ... declare other variables
       LOGICAL :: tlimit, tsteep
@@ -537,7 +533,7 @@
           EXIT DIIS_LOOP
         END IF
 
-        CALL kspotential( 1, .FALSE., tforce, tstress, rhoe, desc, &
+        CALL kspotential( 1, .FALSE., tforce, tstress, rhoe, &
           atoms, bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, fi, vpot, edft, timepre )
 
         s0 = cclock()
@@ -582,7 +578,7 @@
 ! ...       so on).
 
             CALL dforce_all( ispin, c0(:,:,1,ispin), cdesc, fi(:,1,ispin), cgrad(:,:,1,ispin), &
-                             vpot(:,:,:,ispin), eigr, bec )
+                             vpot(:,ispin), eigr, bec )
 
             IF(.NOT.kp%gamma_only) THEN
               DO ik = 1, kp%nkpt
@@ -695,7 +691,6 @@
         USE constants, ONLY: au
         USE cell_base, ONLY: tpiba2
         USE electrons_module, ONLY: eigs, ei, pmss, emass, nb_l, ib_owner, ib_local
-        USE descriptors_module, ONLY: get_local_dims, owner_of, local_index
         USE forces, ONLY: dforce_all
         USE brillouin, ONLY: kpoints, kp
         USE orthogonalize
@@ -714,7 +709,7 @@
         COMPLEX(DP), INTENT(inout) ::  c(:,:,:,:)
         COMPLEX(DP), INTENT(inout) ::  eforce(:,:,:,:)
         TYPE (wave_descriptor), INTENT(in) :: cdesc
-        REAL (DP), INTENT(in) ::  vpot(:,:,:,:), fi(:,:,:)
+        REAL (DP), INTENT(in) ::  vpot(:,:), fi(:,:,:)
         REAL (DP) ::  bec(:,:)
         LOGICAL, INTENT(IN) :: TORTHO
         COMPLEX(DP) :: eigr(:,:)
@@ -749,7 +744,7 @@
 
 ! ...     Calculate | dH / dpsi(j) >
           CALL dforce_all( ispin, c(:,:,1,ispin), cdesc, fi(:,1,ispin), eforce(:,:,1,ispin), &
-                           vpot(:,:,:,ispin), eigr, bec )
+                           vpot(:,ispin), eigr, bec )
 
           DO ik = 1, kp%nkpt
 

CPV/runsd.f90

@@ -28,7 +28,7 @@
 !  -----------------------------------------------------------------------
 !  BEGIN manual
 
-      SUBROUTINE runsd(tortho, tprint, tforce, rhoe, desc, atoms_0, &
+      SUBROUTINE runsd(tortho, tprint, tforce, rhoe, atoms_0, &
                  bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht0, occ, ei, &
                  vpot, doions, edft, maxnstep, sdthr )
 
@@ -47,7 +47,6 @@
       USE atoms_type_module,    ONLY: atoms_type
       USE runcp_module,         ONLY: runcp
       USE phase_factors_module, ONLY: strucf, phfacs
-      USE charge_types,         ONLY: charge_descriptor
       USE control_flags,        ONLY: force_pairing
       use grid_dimensions,      only: nr1, nr2, nr3
       USE reciprocal_vectors,   ONLY: mill_l
@@ -61,9 +60,8 @@
       TYPE (atoms_type), INTENT(INOUT) :: atoms_0
       COMPLEX(DP), INTENT(INOUT) :: c0(:,:,:,:), cm(:,:,:,:), cp(:,:,:,:)
       TYPE (wave_descriptor) :: cdesc
-      REAL(DP) :: rhoe(:,:,:,:)
+      REAL(DP) :: rhoe(:,:)
       COMPLEX(DP) :: sfac(:,:)
-      TYPE (charge_descriptor) :: desc
       COMPLEX(DP) :: eigr(:,:)
       COMPLEX(DP) :: ei1(:,:)
       COMPLEX(DP) :: ei2(:,:)
@@ -75,7 +73,7 @@
       TYPE (dft_energy_type) :: edft
 
       REAL(DP)    :: ei(:,:,:)
-      REAL(DP)    :: vpot(:,:,:,:)
+      REAL(DP)    :: vpot(:,:)
 
       INTEGER   :: maxnstep   !  maximum number of iteration
       REAL(DP) :: sdthr      !  threshold for convergence 
@@ -127,7 +125,7 @@
 
         s1 = cclock()
 
-        CALL kspotential( 1, ttprint, ttforce, ttstress, rhoe, desc, atoms_0, &
+        CALL kspotential( 1, ttprint, ttforce, ttstress, rhoe, atoms_0, &
                           bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0,  &
                           occ, vpot, edft, timepre )
 
@@ -171,7 +169,7 @@
 
       IF( tforce ) THEN
         atoms_0%for = 0.0d0
-        CALL kspotential( 1, ttprint, tforce, ttstress, rhoe, desc, &
+        CALL kspotential( 1, ttprint, tforce, ttstress, rhoe, &
           atoms_0, bec, becdr, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, vpot, edft, timepre )
         IF(ionode ) THEN
           WRITE( stdout,fmt="(12X,'runsd: fion and edft calculated = ',F14.6)") edft%etot

CPV/smcp.f90

@@ -93,6 +93,8 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
   USE mp_global,                ONLY : mp_global_start
   USE mp,                       ONLY : mp_sum
   USE fft_base,                 ONLY : dfftp
+  USE orthogonalize,            ONLY : ortho
+  USE orthogonalize_base,       ONLY : updatc, calphi
   !
 #if ! defined __NOSMD
   !
@@ -612,8 +614,8 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
         !     imposing the orthogonality
         !     ==========================================================
         !
-        CALL calphi( rep_el(sm_k)%cm, ngw, ema0bg,rep_el(sm_k)%bec, nkb, &
+        CALL calphi( rep_el(sm_k)%cm, ngw, rep_el(sm_k)%bec, nkb, &
-             & vkb,rep_el(sm_k)%phi, nbsp )
+             & vkb,rep_el(sm_k)%phi, nbsp, ema0bg )
         !
         !
         IF(ionode) WRITE( sm_file,*) ' out from calphi'
@@ -621,7 +623,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
         !
         IF(tortho) THEN
            CALL ortho  (eigr,rep_el(sm_k)%c0,rep_el(sm_k)%phi,rep_el(sm_k)%lambda, &
-                &                   bigr,iter,ccc(sm_k),ortho_eps,ortho_max,delt,bephi,becp)
+                &                   bigr,iter,ccc(sm_k),bephi,becp)
         ELSE
            CALL gram( vkb, rep_el(sm_k)%bec, nkb, rep_el(sm_k)%c0, ngw, nbsp )
            !
@@ -643,8 +645,9 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
         ENDIF
         !
         IF(tortho) THEN
-           CALL updatc(ccc(sm_k),rep_el(sm_k)%lambda,rep_el(sm_k)%phi, &
+           CALL updatc( ccc(sm_k), nbsp, rep_el(sm_k)%lambda, SIZE( rep_el(sm_k)%lambda, 1 ), &
-                & bephi,becp,rep_el(sm_k)%bec,rep_el(sm_k)%c0)
+                        rep_el(sm_k)%phi, SIZE( rep_el(sm_k)%phi, 1 ), bephi, SIZE(bephi,1), &
+                        becp,rep_el(sm_k)%bec,rep_el(sm_k)%c0)
            !
            IF(ionode) WRITE( sm_file,*) ' out from updatc'
         ENDIF
@@ -959,7 +962,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
         !     calphi calculates phi
         !     the electron mass rises with g**2
         !
-        CALL calphi( rep_el(sm_k)%c0, ngw, ema0bg, rep_el(sm_k)%bec, nkb, vkb, rep_el(sm_k)%phi, nbsp )
+        CALL calphi( rep_el(sm_k)%c0, ngw, rep_el(sm_k)%bec, nkb, vkb, rep_el(sm_k)%phi, nbsp, ema0bg )
         !
         !     begin try and error loop (only one step!)
         !
@@ -1171,7 +1174,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
         IF(tortho) THEN
            CALL ortho                                                     &
                 &         (eigr,rep_el(sm_k)%cm,rep_el(sm_k)%phi,rep_el(sm_k)%lambda, &
-                & bigr,iter,ccc(sm_k),ortho_eps,ortho_max,delt,bephi,becp)
+                & bigr,iter,ccc(sm_k),bephi,becp)
         ELSE
            CALL gram( vkb, rep_el(sm_k)%bec, nkb, rep_el(sm_k)%cm, ngw, nbsp )
            IF(iprsta.GT.4) CALL dotcsc(eigr,rep_el(sm_k)%cm)
@@ -1183,8 +1186,10 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
         !
         IF(iprsta.GE.3) CALL print_lambda( rep_el(sm_k)%lambda, nbsp, 9, 1.0d0 )
         !
-        IF(tortho) CALL updatc(ccc(sm_k),rep_el(sm_k)%lambda,rep_el(sm_k)%phi,bephi, &
+        IF(tortho) &
-             & becp,rep_el(sm_k)%bec,rep_el(sm_k)%cm)
+           CALL updatc( ccc(sm_k), nbsp, rep_el(sm_k)%lambda, SIZE(rep_el(sm_k)%lambda,1), &
+                        rep_el(sm_k)%phi, SIZE(rep_el(sm_k)%phi,1), bephi, SIZE(bephi,1), &
+                        becp, rep_el(sm_k)%bec, rep_el(sm_k)%cm )
         !
         CALL calbec (nvb+1,nsp,eigr,rep_el(sm_k)%cm,rep_el(sm_k)%bec)
         IF (tpre) CALL caldbec(ngw,nkb,nbsp,1,nsp,eigr,rep_el(sm_k)%cm,dbec,.true.)

CPV/stress.f90

@@ -85,7 +85,7 @@
 
 ! ... declare subroutine arguments
       REAL(DP) :: pail(:,:), desr(:), strvxc
-      REAL(DP) :: grho(:,:,:,:,:), v2xc(:,:,:,:,:)
+      REAL(DP) :: grho(:,:,:), v2xc(:,:,:)
       REAL(DP) :: bec(:,:)
       COMPLEX(DP) :: rhoeg(:,:), vxc(:,:)
       COMPLEX(DP), INTENT(IN) :: sfac(:,:)

CPV/turbo.f90

@@ -15,7 +15,7 @@
 
         LOGICAL :: TTURBO
         INTEGER :: NTURBO
-        COMPLEX(DP), ALLOCATABLE :: turbo_states(:,:,:,:) 
+        COMPLEX(DP), ALLOCATABLE :: turbo_states(:,:) 
 
         PUBLIC :: tturbo, nturbo, turbo_states, turbo_init, allocate_turbo
         PUBLIC :: deallocate_turbo
@@ -35,19 +35,19 @@
           RETURN
         END SUBROUTINE turbo_init
 
-        SUBROUTINE allocate_turbo( nr1, nr2, nr3 )
+        SUBROUTINE allocate_turbo( nnr )
           USE io_global, ONLY: ionode
           USE io_global, ONLY: stdout
 
           USE mp, ONLY: mp_sum
-          INTEGER :: nr1,nr2,nr3
+          INTEGER :: nnr
           INTEGER :: ierr
           IF( ionode ) THEN
             WRITE( stdout,fmt='(/,3X,"TURBO: allocating ",I10," bytes ")') &
-              16*nr1*nr2*nr3*nturbo
+              16*nnr*nturbo
           END IF
           IF( .NOT. ALLOCATED( turbo_states ) ) THEN
-            ALLOCATE( turbo_states( nr1, nr2, nr3, nturbo ), STAT = ierr)
+            ALLOCATE( turbo_states( nnr, nturbo ), STAT = ierr)
             CALL mp_sum(ierr)
             IF( ierr /= 0 ) THEN 
               IF( ionode ) THEN

CPV/wave.f90

@@ -7,27 +7,11 @@
 !
 #include "f_defs.h"
 
-!  ----------------------------------------------
-!  BEGIN manual
 
 !=----------------------------------------------------------------------------=!
    MODULE wave_functions
 !=----------------------------------------------------------------------------=!
 
-!  (describe briefly what this module does...)
-!  ----------------------------------------------
-!  routines in this module:
-!  REAL(DP) FUNCTION dft_kinetic_energy(c,hg,f,nb)
-!  REAL(DP) FUNCTION cp_kinetic_energy(cp,cm,pmss,emass,delt)
-!  SUBROUTINE update_wave_functions(cm,c0,cp)
-!  SUBROUTINE crot_gamma (c0,lambda,eig)
-!  SUBROUTINE crot_kp (ik,c0,lambda,eig)
-!  SUBROUTINE proj_gamma(a,b,lambda)
-!  SUBROUTINE proj_kp(ik,a,b,lambda)
-!  ----------------------------------------------
-!  END manual
-
-
 ! ...   include modules
         USE kinds
 
@@ -47,202 +31,13 @@
             MODULE PROCEDURE fixwave_s, fixwave_v, fixwave_m
           END INTERFACE
 
-          PUBLIC :: dft_kinetic_energy, cp_kinetic_energy
+          PUBLIC :: cp_kinetic_energy
           PUBLIC :: update_wave_functions, wave_rand_init
 
-!  end of module-scope declarations
-!  ----------------------------------------------
-
 !=----------------------------------------------------------------------------=!
       CONTAINS
 !=----------------------------------------------------------------------------=!
 
-!  subroutines
-!  ----------------------------------------------
-!  ----------------------------------------------
-
-
-    REAL(DP) FUNCTION dft_kinetic_energy(c0, cdesc, f, xmkin)
-
-!  This function compute the Total Quanto-Mechanical Kinetic Energy of the Kohn-Sham
-!  wave function
-!  ----------------------------------------------
-
-        USE cell_base, ONLY: tpiba2
-        USE brillouin, ONLY: kpoints, kp
-        USE wave_types, ONLY: wave_descriptor
-        USE electrons_module, ONLY: pmss
-        USE control_flags, ONLY: force_pairing, gamma_only
-        USE reciprocal_space_mesh, ONLY: gkcutz_l
-        USE reciprocal_vectors,    ONLY: ggp
-
-        IMPLICIT NONE
-
-        COMPLEX(DP), INTENT(IN) :: c0(:,:,:,:)       !  wave functions coefficients
-        TYPE (wave_descriptor), INTENT(IN) :: cdesc   !  descriptor of c0
-        REAL(DP), INTENT(IN) :: f(:,:,:)             !  occupation numbers
-        REAL(DP), OPTIONAL, INTENT(INOUT) :: xmkin
-
-        INTEGER    :: ib, ik, ispin, ispin_wfc
-        REAL(DP)  :: xkin, fact, xkins, xmkins
-
-! ... end of declarations
-!  ----------------------------------------------
-
-
-      xkin = 0.d0
-
-      DO ispin = 1, cdesc%nspin
-
-        ispin_wfc = ispin
-        IF( force_pairing ) ispin_wfc = 1
-         
-        DO ik = 1, cdesc%nkl
-
-          fact = kp%weight(ik)
-          IF( cdesc%gamma ) THEN
-            fact = fact * 2.d0
-          END IF
-
-          IF( cdesc%gamma ) THEN
-             xkins = dft_kinetic_energy_s( ispin, c0(:,:,ik,ispin_wfc), cdesc, ggp, f(:,ik,ispin) )
-             IF( PRESENT( xmkin ) ) THEN
-               xmkins = dft_weighted_kinene( ispin, c0(:,:,ik,ispin_wfc), cdesc, ggp, f(:,ik,ispin) )
-             END IF
-          ELSE
-             xkins = dft_kinetic_energy_s( ispin, c0(:,:,ik,ispin_wfc), cdesc, gkcutz_l(:,ik), f(:,ik,ispin) )
-             IF( PRESENT( xmkin ) ) THEN
-               xmkins =  dft_weighted_kinene( ispin, c0(:,:,ik,ispin_wfc), cdesc, gkcutz_l(:,ik), f(:,ik,ispin) )
-             END IF
-          ENDIF
-
-          xkin  = xkin + fact * xkins
-
-          IF( PRESENT( xmkin ) ) THEN
-            xmkin = xmkin + fact * xmkins
-          END IF
-
-        END DO
-
-      END DO
-
-      dft_kinetic_energy = xkin
-
-      RETURN
-    END FUNCTION dft_kinetic_energy
-
-!=----------------------------------------------------------------------------=!
-
-
-      REAL(DP) FUNCTION dft_weighted_kinene( ispin, c, cdesc, g2, fi)
-
-!  (describe briefly what this routine does...)
-!  ----------------------------------------------
-
-        USE wave_types, ONLY: wave_descriptor
-        USE electrons_module, ONLY: pmss
-
-        COMPLEX(DP), INTENT(IN) :: c(:,:)
-        INTEGER, INTENT( IN ) :: ispin
-        TYPE (wave_descriptor), INTENT(IN) :: cdesc
-        REAL (DP), INTENT(IN) :: fi(:), g2(:)
-        INTEGER    ib, ig
-        REAL(DP)  skm, xmkin
-! ... end of declarations
-
-        xmkin = 0.0d0
-
-        IF( cdesc%nbl( ispin ) > SIZE( c, 2 ) .OR. &
-            cdesc%nbl( ispin ) > SIZE( fi )        ) &
-          CALL errore( ' dft_weighted_kinene ', ' wrong sizes ', 1 )
-        IF( cdesc%ngwl > SIZE( c, 1 ) .OR. &
-            cdesc%ngwl > SIZE( g2 )     .OR. &
-            cdesc%ngwl > SIZE( pmss )      ) &
-          CALL errore( ' dft_weighted_kinene ', ' wrong sizes ', 2 )
-
-        IF( cdesc%gamma .AND. cdesc%gzero ) THEN
-
-          DO ib = 1, cdesc%nbl( ispin )
-            skm = 0.d0
-            DO ig = 2, cdesc%ngwl
-              skm  = skm + g2(ig) * DBLE( CONJG(c(ig,ib)) * c(ig,ib) ) * pmss(ig)
-            END DO
-            skm = skm + g2(1) * DBLE( c(1,ib) )**2 * pmss(1) / 2.0d0
-            xmkin = xmkin + fi(ib) * skm * 0.5d0
-          END DO
-
-        ELSE
-
-          DO ib = 1, cdesc%nbl( ispin )
-            skm = 0.d0
-            DO ig = 1, cdesc%ngwl
-              skm  = skm + g2(ig) * DBLE( CONJG( c( ig, ib ) ) * c( ig, ib ) ) * pmss(ig)
-            END DO
-            xmkin = xmkin + fi(ib) * skm * 0.5d0
-          END DO
-
-        END IF
-
-        dft_weighted_kinene = xmkin
-
-        RETURN
-      END FUNCTION dft_weighted_kinene
-
-!=----------------------------------------------------------------------------=!
-
-      REAL(DP) FUNCTION dft_kinetic_energy_s( ispin, c, cdesc, g2, fi)
-
-!  (describe briefly what this routine does...)
-!  ----------------------------------------------
-
-        USE wave_types, ONLY: wave_descriptor
-        COMPLEX(DP), INTENT(IN) :: c(:,:)
-        INTEGER, INTENT( IN ) :: ispin
-        TYPE (wave_descriptor), INTENT(IN) :: cdesc
-        REAL (DP),  INTENT(IN) :: fi(:), g2(:)
-        INTEGER    ib, ig, igs
-        REAL(DP)  sk1, xkin
-! ... end of declarations
-
-        xkin = 0.0d0
-
-        IF( cdesc%nbl( ispin ) > SIZE( c, 2 ) .OR. &
-            cdesc%nbl( ispin ) > SIZE( fi )        ) &
-          CALL errore( ' dft_total_charge ', ' wrong sizes ', 1 )
-        IF( cdesc%ngwl > SIZE( c, 1 ) .OR. &
-            cdesc%ngwl > SIZE( g2 )      ) &
-          CALL errore( ' dft_total_charge ', ' wrong sizes ', 2 )
-
-        IF( cdesc%gamma .AND. cdesc%gzero ) THEN
-
-          DO ib = 1, cdesc%nbl( ispin )
-            sk1 = 0.d0
-            DO ig = 2, cdesc%ngwl
-              sk1 = sk1 + g2(ig) * DBLE( CONJG( c(ig,ib) ) * c(ig,ib) )
-            END DO
-            sk1  = sk1 + g2(1) * DBLE( c(1,ib) )**2 / 2.0d0
-            xkin = xkin + fi(ib) * sk1 * 0.5d0
-          END DO
-
-        ELSE
-
-          DO ib = 1, cdesc%nbl( ispin )
-            sk1 = 0.d0
-            DO ig = 1, cdesc%ngwl
-              sk1 = sk1 + g2(ig) * DBLE( CONJG( c(ig,ib) ) * c(ig,ib) )
-            END DO
-            xkin = xkin + fi(ib) * sk1 * 0.5d0
-          END DO
-
-        END IF
-
-        dft_kinetic_energy_s = xkin
-
-        RETURN
-      END FUNCTION dft_kinetic_energy_s
-
-
-!=----------------------------------------------------------------------------=!
 
    SUBROUTINE fixwave_s ( ispin, c, cdesc, kmask )
 

CPV/wf.f90

@@ -3569,7 +3569,7 @@ SUBROUTINE dforce_field( bec, deeq, betae, i, c, ca, df, da, v, v1 )
   USE cvan,                   ONLY : ish
   USE smooth_grid_dimensions, ONLY : nr1s, nr2s, nr3s, &
                                      nr1sx, nr2sx, nr3sx, nnrsx
-  USE electrons_base,         ONLY : nbspx, nbsp, nspin, f, fspin
+  USE electrons_base,         ONLY : nbspx, nbsp, nspin, f, ispin
   USE constants,              ONLY : pi, fpi
   USE ions_base,              ONLY : nsp, na, nat
   USE gvecw,                  ONLY : ggp
@@ -3626,12 +3626,12 @@ SUBROUTINE dforce_field( bec, deeq, betae, i, c, ca, df, da, v, v1 )
      !
   ENDIF
   ! 
-  iss1=fspin(i)
+  iss1=ispin(i)
   !
   ! the following avoids a potential out-of-bounds error
   !
   IF (i.NE.nbsp) THEN
-     iss2=fspin(i+1)
+     iss2=ispin(i+1)
   ELSE
      iss2=iss1
   END IF
@@ -3874,7 +3874,7 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
   USE cell_base,              ONLY : omega
   USE smooth_grid_dimensions, ONLY : nr1s, nr2s, nr3s, &
                                      nr1sx, nr2sx, nr3sx, nnrsx
-  USE electrons_base,         ONLY : nbspx, nbsp, nspin, f, fspin
+  USE electrons_base,         ONLY : nbspx, nbsp, nspin, f, ispin
   USE constants,              ONLY : pi, fpi
   USE wannier_base,           ONLY : iwf
   USE dener,                  ONLY : dekin, denl
@@ -3914,7 +3914,7 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
   !     ==================================================================
   !     calculation of kinetic energy ekin
   !     ==================================================================
-  ekin=enkin(c)
+  ekin=enkin(c,ngw,f,nbsp)
   IF(tpre) CALL denkin(c,dekin)
   !
   !     ==================================================================
@@ -4018,11 +4018,11 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
 #else
      IF(tbuff) WRITE(21,iostat=ios) psis
 #endif
-     !            iss1=fspin(i)
+     !            iss1=ispin(i)
      iss1=1
      sa1=f(i)/omega
      !            if (i.ne.nbsp) then
-     !              iss2=fspin(i+1)
+     !              iss2=ispin(i+1)
      !              sa2=f(i+1)/omega
      !            else
      iss2=iss1  ! carlo

Modules/descriptors.f90

@@ -12,84 +12,16 @@
         IMPLICIT NONE
         SAVE
 
-        INTERFACE desc_init
+        INTEGER  ldim_block, ldim_cyclic, ldim_block_cyclic
-          MODULE PROCEDURE desc_init_1d, desc_init_2d, desc_init_3d
+        INTEGER  lind_block, lind_cyclic, lind_block_cyclic
-        END INTERFACE
+        EXTERNAL ldim_block, ldim_cyclic, ldim_block_cyclic
-        INTERFACE global_index
+        EXTERNAL lind_block, lind_cyclic, lind_block_cyclic
-          MODULE PROCEDURE globalindex_desc, globalindex_shape
-        END INTERFACE
-        INTERFACE local_index
-          MODULE PROCEDURE localindex_desc, localindex_shape
-        END INTERFACE
-        INTERFACE local_dimension
-          MODULE PROCEDURE localdim_desc, localdim_shape
-        END INTERFACE
-        INTERFACE owner_of
-          MODULE PROCEDURE ownerof_desc, ownerof_shape
-        END INTERFACE
-        INTERFACE get_local_dims
-          MODULE PROCEDURE desc_ldims
-        END INTERFACE
-        INTERFACE get_global_dims
-          MODULE PROCEDURE desc_gdims
-        END INTERFACE
-
-        INTEGER NUMROC
-        EXTERNAL NUMROC
        
         CONTAINS
 
 
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-!
-          SUBROUTINE desc_init_1d(desc, matrix_type, rows, &
-            row_block, row_src_pe, grid, row_shape)
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-              TYPE (descriptor) :: desc
-              INTEGER, INTENT(IN) :: matrix_type    
-              TYPE (processors_grid), INTENT(IN) :: grid 
-              INTEGER, INTENT(IN) :: rows   
-              INTEGER, INTENT(IN) :: row_block
-              INTEGER, INTENT(IN) :: row_src_pe   
-              INTEGER, INTENT(IN) :: row_shape   
-
-              desc%matrix_type = matrix_type
-              desc%grid = grid
-
-              CALL desc_init_x(desc%nx, desc%xshape, desc%nxl, &
-                desc%nxblk, desc%ixl, desc%ipexs, rows, row_shape, &
-                row_block, row_src_pe, grid%mex, grid%npx)
-
-              desc%ny = 1
-              desc%nyl = 1
-              desc%nyblk = 1
-              desc%ipeys = 0
-              desc%yshape = REPLICATED_DATA_SHAPE
-
-              desc%nz = 1
-              desc%nzl = 1
-              desc%nzblk = 1
-              desc%ipezs = 0
-              desc%zshape = REPLICATED_DATA_SHAPE
-
-              desc%ldx = 1
-              desc%ldy = 1
-            RETURN
-          END SUBROUTINE desc_init_1d
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
           SUBROUTINE desc_init_blacs(desc, matrix_type, rows, columns, &
             row_block, column_block, row_src_pe, column_src_pe, grid, local_ld)
-!
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
 
               TYPE (descriptor) :: desc
               INTEGER, INTENT(IN) :: matrix_type    
@@ -105,201 +37,73 @@
               desc%matrix_type = matrix_type
               desc%grid = grid
 
-              CALL desc_init_x(desc%nx, desc%xshape, desc%nxl, &
+              CALL desc_init_x(desc%nx, desc%xdist, desc%nxl, &
                 desc%nxblk, desc%ixl, desc%ipexs, rows, &
-                BLOCK_CYCLIC_SHAPE, row_block, row_src_pe, grid%mex, &
+                BLOCK_CYCLIC_DIST, row_block, row_src_pe, grid%mex, &
                 grid%npx)
-              CALL desc_init_x(desc%ny, desc%yshape, &
+              CALL desc_init_x(desc%ny, desc%ydist, &
                 desc%nyl, desc%nyblk, desc%iyl, &
-                desc%ipeys, columns, BLOCK_CYCLIC_SHAPE, column_block, &
+                desc%ipeys, columns, BLOCK_CYCLIC_DIST, column_block, &
                 column_src_pe, grid%mey, grid%npy)
 
               desc%nz = 1
               desc%nzl = 1
               desc%nzblk = 1
               desc%ipezs = 0
-              desc%zshape = REPLICATED_DATA_SHAPE
+              desc%zdist = REPLICATED_DATA_DIST
 
               IF(PRESENT(local_ld)) THEN
                 desc%ldx = local_ld
               ELSE
-                desc%ldx = localdim_shape( rows, row_block, grid%mex, &
+                desc%ldx = ldim_block_cyclic( rows, row_block, grid%npx, grid%mex )
-                  row_src_pe, grid%npx, desc%xshape)
               END IF
               desc%ldy = 1
 
             RETURN
           END SUBROUTINE desc_init_blacs
 
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
 
-          SUBROUTINE desc_init_2d(desc, matrix_type, rows, columns, &
+          SUBROUTINE desc_init_x(desc_nxs, desc_nx_dist, desc_local_nxs, &
-            row_block, column_block, row_src_pe, &
+            desc_nx_block, desc_ix, desc_nx_src_pe, nxs, nx_dist, nx_block, &
-            column_src_pe, grid, row_shape, column_shape, local_ld)
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-              TYPE (descriptor) :: desc
-              INTEGER, INTENT(IN) :: matrix_type    
-              TYPE (processors_grid), INTENT(IN) :: grid 
-              INTEGER, INTENT(IN) :: rows   
-              INTEGER, INTENT(IN) :: columns
-              INTEGER, INTENT(IN) :: row_block
-              INTEGER, INTENT(IN) :: column_block
-              INTEGER, INTENT(IN) :: row_src_pe   
-              INTEGER, INTENT(IN) :: column_src_pe
-              INTEGER, INTENT(IN) :: row_shape   
-              INTEGER, INTENT(IN) :: column_shape
-              INTEGER, INTENT(IN), OPTIONAL :: local_ld
- 
-              LOGICAL :: debug = .FALSE.
-
-              desc%matrix_type = matrix_type
-              desc%grid = grid
-
-              CALL desc_init_x(desc%nx, desc%xshape, desc%nxl, &
-                desc%nxblk, desc%ixl, desc%ipexs, rows, row_shape, &
-                row_block, row_src_pe, grid%mex, grid%npx)
-              CALL desc_init_x(desc%ny, desc%yshape, &
-                desc%nyl, desc%nyblk, desc%iyl, &
-                desc%ipeys, columns, column_shape, column_block, &
-                column_src_pe, grid%mey, grid%npy)
-
-              IF( debug ) THEN
-                WRITE( stdout,fmt="(' desc%nx      = ', I6 )") desc%nx
-                WRITE( stdout,fmt="(' desc%xshape  = ', I6 )") desc%xshape
-                WRITE( stdout,fmt="(' desc%nxl     = ', I6 )") desc%nxl
-                WRITE( stdout,fmt="(' desc%nxblk   = ', I6 )") desc%nxblk
-                WRITE( stdout,fmt="(' desc%ixl     = ', I6 )") desc%ixl
-                WRITE( stdout,fmt="(' desc%ipexs   = ', I6 )") desc%ipexs
-
-                WRITE( stdout,fmt="(' desc%ny      = ', I6 )") desc%ny
-                WRITE( stdout,fmt="(' desc%yshape  = ', I6 )") desc%yshape
-                WRITE( stdout,fmt="(' desc%nyl     = ', I6 )") desc%nyl
-                WRITE( stdout,fmt="(' desc%nyblk   = ', I6 )") desc%nyblk
-                WRITE( stdout,fmt="(' desc%iyl     = ', I6 )") desc%iyl
-                WRITE( stdout,fmt="(' desc%ipeys   = ', I6 )") desc%ipeys
-              END IF
-
-              desc%nz = 1
-              desc%nzl = 1
-              desc%nzblk = 1
-              desc%ipezs = 0
-              desc%zshape = REPLICATED_DATA_SHAPE
-
-              IF(PRESENT(local_ld)) THEN
-                desc%ldx = local_ld
-              ELSE
-                desc%ldx = localdim_shape( rows, row_block, grid%mex, &
-                  row_src_pe, grid%npx, desc%xshape)
-              END IF
-              desc%ldy = 1
-            RETURN
-          END SUBROUTINE desc_init_2d
-
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          SUBROUTINE desc_init_3d(desc, matrix_type, rows, columns, &
-            planes, row_block, column_block, plane_block, row_src_pe, &
-            column_src_pe, plane_src_pe, grid, row_shape, column_shape, &
-            plane_shape, local_ld, local_sub_ld)
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-              TYPE (descriptor) :: desc
-              INTEGER, INTENT(IN) :: matrix_type    
-              TYPE (processors_grid), INTENT(IN) :: grid 
-              INTEGER, INTENT(IN) :: rows   
-              INTEGER, INTENT(IN) :: columns
-              INTEGER, INTENT(IN) :: planes
-              INTEGER, INTENT(IN) :: row_block
-              INTEGER, INTENT(IN) :: column_block
-              INTEGER, INTENT(IN) :: plane_block
-              INTEGER, INTENT(IN) :: row_src_pe   
-              INTEGER, INTENT(IN) :: column_src_pe
-              INTEGER, INTENT(IN) :: plane_src_pe
-              INTEGER, INTENT(IN) :: row_shape   
-              INTEGER, INTENT(IN) :: column_shape
-              INTEGER, INTENT(IN) :: plane_shape
-              INTEGER, INTENT(IN), OPTIONAL :: local_ld
-              INTEGER, INTENT(IN), OPTIONAL :: local_sub_ld
-
-              desc%matrix_type = matrix_type
-              desc%grid = grid
-
-              CALL desc_init_x(desc%nx, desc%xshape, desc%nxl, &
-                desc%nxblk, desc%ixl, desc%ipexs, rows, row_shape, &
-                row_block, row_src_pe, grid%mex, grid%npx)
-              CALL desc_init_x(desc%ny, desc%yshape, &
-                desc%nyl, desc%nyblk, desc%iyl, &
-                desc%ipeys, columns, column_shape, column_block, &
-                column_src_pe, grid%mey, grid%npy)
-              CALL desc_init_x(desc%nz, desc%zshape, &
-                desc%nzl, desc%nzblk, desc%izl, &
-                desc%ipezs, planes, plane_shape, plane_block, &
-                plane_src_pe, grid%mez, grid%npz)
-
-              IF(PRESENT(local_ld)) THEN
-                desc%ldx = local_ld
-              ELSE
-                desc%ldx = localdim_shape( rows, row_block, grid%mex, &
-                  row_src_pe, grid%npx, desc%xshape)
-              END IF
-              IF(PRESENT(local_sub_ld)) THEN
-                desc%ldy = local_sub_ld
-              ELSE
-                desc%ldy =  localdim_shape( columns, column_block, &
-                  grid%mey, column_src_pe, grid%npy, &
-                  desc%yshape)
-              END IF
-            RETURN
-          END SUBROUTINE desc_init_3d
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          SUBROUTINE desc_init_x(desc_nxs, desc_nx_shape, desc_local_nxs, &
-            desc_nx_block, desc_ix, desc_nx_src_pe, nxs, nx_shape, nx_block, &
             nx_src_pe, mype, npes)
 !
-!  END manual
-!=----------------------------------------------------------------------------=!
 
               IMPLICIT NONE
               INTEGER, INTENT(OUT) ::  desc_nxs
-              INTEGER, INTENT(OUT) ::  desc_nx_shape
+              INTEGER, INTENT(OUT) ::  desc_nx_dist
               INTEGER, INTENT(OUT) ::  desc_local_nxs
               INTEGER, INTENT(OUT) ::  desc_nx_block
               INTEGER, INTENT(OUT) ::  desc_ix
               INTEGER, INTENT(OUT) ::  desc_nx_src_pe
               INTEGER, INTENT(IN)  ::  nxs
-              INTEGER, INTENT(IN)  ::  nx_shape
+              INTEGER, INTENT(IN)  ::  nx_dist
               INTEGER, INTENT(IN)  ::  nx_block
               INTEGER, INTENT(IN)  ::  nx_src_pe
               INTEGER, INTENT(IN)  ::  mype
               INTEGER, INTENT(IN)  ::  npes
                 
               desc_nxs      = nxs
-              desc_nx_shape  = nx_shape
+              desc_nx_dist  = nx_dist
-              desc_local_nxs = localdim_shape( nxs, nx_block, mype, nx_src_pe, npes, nx_shape)
-              desc_ix        = localindex_shape( 1, nxs, nx_block, mype, npes, nx_shape)
 
-              SELECT CASE (nx_shape)
+              SELECT CASE (nx_dist)
-              CASE ( BLOCK_CYCLIC_SHAPE )
+              CASE ( BLOCK_CYCLIC_DIST )
+                desc_local_nxs = ldim_block_cyclic( nxs, nx_block, npes, mype )
+                desc_ix        = lind_block_cyclic( 1, nxs, nx_block, npes, mype)
                 desc_nx_block  = nx_block
                 desc_nx_src_pe = nx_src_pe
-              CASE ( BLOCK_PARTITION_SHAPE )
+              CASE ( BLOCK_PARTITION_DIST )
+                desc_local_nxs = ldim_block( nxs, npes, mype )
+                desc_ix        = lind_block( 1, nxs, npes, mype)
                 desc_nx_block  = desc_local_nxs 
                 desc_nx_src_pe = 0             
-              CASE ( CYCLIC_SHAPE )
+              CASE ( CYCLIC_DIST )
+                desc_local_nxs = ldim_cyclic( nxs, npes, mype )
+                desc_ix        = lind_cyclic( 1, nxs, npes, mype)
                 desc_nx_block  = 1 
                 desc_nx_src_pe = 0             
-              CASE ( REPLICATED_DATA_SHAPE )
+              CASE ( REPLICATED_DATA_DIST )
+                desc_local_nxs = nxs
+                desc_ix        = 1
                 desc_nx_block  = nxs
                 desc_nx_src_pe = mype
               END SELECT
@@ -329,324 +133,5 @@
             RETURN
           END SUBROUTINE pblas_descriptor
 
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          INTEGER FUNCTION globalindex_shape( lind, n, nb, me, isrc, np, pshape )
-
-!  This function computes the global index of a distributed array entry
-!  pointed to by the local index lind of the process indicated by me.
-!  lind      local index of the distributed matrix entry.
-!  N         is the size of the global array.
-!  NB        size of the blocks the distributed matrix is split into.
-!  me        The coordinate of the process whose local array row or
-!            column is to be determined.
-!  isrc      The coordinate of the process that possesses the first
-!            row/column of the distributed matrix.
-!  np        The total number processes over which the distributed
-!            matrix is distributed.
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-
-            INTEGER, INTENT(IN) :: lind, n, nb, me, isrc, np, pshape
-            INTEGER r, q
-
-            IF( pshape .EQ. BLOCK_PARTITION_SHAPE ) THEN
-
-              q = INT(n/np)
-              r = MOD(n,np)
-              IF( me < r ) THEN
-                GLOBALINDEX_SHAPE = (Q+1)*me + lind
-              ELSE
-                GLOBALINDEX_SHAPE = Q*me + R + lind
-              END IF
-
-            ELSE IF ( pshape .EQ. BLOCK_CYCLIC_SHAPE ) THEN
-
-              GLOBALINDEX_SHAPE = np*NB*((lind-1)/NB) + &
-                MOD(lind-1,NB) + MOD(np+me-isrc, np)*NB + 1
-
-            ELSE IF ( pshape .EQ. CYCLIC_SHAPE ) THEN
-
-                GLOBALINDEX_SHAPE = (lind-1) * np + me + 1
-
-            ELSE
-
-              GLOBALINDEX_SHAPE = lind
-
-            END IF
-            RETURN
-          END FUNCTION globalindex_shape
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          INTEGER FUNCTION globalindex_desc( lind, desc, what )
-
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-            INTEGER, INTENT(IN) :: lind
-            TYPE (descriptor) :: desc
-            CHARACTER(LEN=*) :: what
-            INTEGER N, nb, src_pe, my_pe, np, pshape
-            IF ( what(1:1) .EQ. 'R' .OR. what(1:1) .EQ. 'r' ) THEN
-              NB     = desc%nxblk;   N      = desc%nx; 
-              np     = desc%grid%npx; src_pe = desc%ipexs; 
-              my_pe  = desc%grid%mex; pshape = desc%xshape
-            ELSE IF ( what(1:1) .EQ. 'C' .OR. what(1:1) .EQ. 'c' ) THEN
-              NB     = desc%nyblk;   N      = desc%ny; 
-              np     = desc%grid%npy; src_pe = desc%ipeys; 
-              my_pe  = desc%grid%mey; pshape = desc%yshape
-            ELSE IF ( what(1:1) .EQ. 'P' .OR. what(1:1) .EQ. 'p' ) THEN
-              NB     = desc%nzblk;   N      = desc%nz; 
-              np     = desc%grid%npz; src_pe = desc%ipezs; 
-              my_pe  = desc%grid%mez; pshape = desc%zshape
-            END IF
-            globalindex_desc = globalindex_shape(lind, n, nb, my_pe, src_pe, np, pshape )
-            RETURN
-          END FUNCTION globalindex_desc
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          INTEGER FUNCTION localdim_shape( n, nb, me, isrc, np, pshape)
-            
-!  N     = Global dimension of the array
-!  NB    = Size of the blocks ( meaningful only for BLOCK_CYCLIC_SHAPE )
-!  me = Index of the callig processor
-!  isrc = Index of the processor owning the first element of the array
-!  np = Number of processors among which the array is subdivided
-!  pshape  = Shape of the distributed data
-!
-!  This function return the number of array elements owned
-!  by the callig processor
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-           IMPLICIT NONE
-           INTEGER, INTENT(IN) :: n, nb, me, isrc, np, pshape
-
-           IF( pshape .EQ. BLOCK_PARTITION_SHAPE ) THEN
-
-             LOCALDIM_SHAPE = INT(N/np)
-             IF( me < MOD(N,np) ) LOCALDIM_SHAPE = LOCALDIM_SHAPE + 1
-
-           ELSE IF( pshape .EQ. BLOCK_CYCLIC_SHAPE ) THEN
-
-             LOCALDIM_SHAPE = NUMROC( N, NB, me, isrc, np )
-
-           ELSE IF( pshape .EQ. CYCLIC_SHAPE ) THEN
-
-             LOCALDIM_SHAPE = INT(N/np)
-             IF( me < MOD(N,np) ) LOCALDIM_SHAPE = LOCALDIM_SHAPE + 1
-
-           ELSE
-
-             LOCALDIM_SHAPE = n
-
-           END IF
-           RETURN
-         END FUNCTION localdim_shape
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          INTEGER FUNCTION localdim_desc( desc, what )
-
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-            TYPE (descriptor) :: desc
-            CHARACTER(LEN=*) :: what
-            INTEGER n, nb, src_pe, my_pe, np, pshape
-            IF ( what(1:1) .EQ. 'R' .OR. what(1:1) .EQ. 'r' ) THEN
-              NB     = desc%nxblk;   N      = desc%nx; 
-              np     = desc%grid%npx; src_pe = desc%ipexs; 
-              my_pe  = desc%grid%mex; pshape = desc%xshape
-            ELSE IF ( what(1:1) .EQ. 'C' .OR. what(1:1) .EQ. 'c' ) THEN
-              NB     = desc%nyblk;   N      = desc%ny; 
-              np     = desc%grid%npy; src_pe = desc%ipeys; 
-              my_pe  = desc%grid%mey; pshape = desc%yshape
-            ELSE IF ( what(1:1) .EQ. 'P' .OR. what(1:1) .EQ. 'p' ) THEN
-              NB     = desc%nzblk;   N      = desc%nz; 
-              np     = desc%grid%npz; src_pe = desc%ipezs; 
-              my_pe  = desc%grid%mez; pshape = desc%zshape
-            END IF
-            localdim_desc = localdim_shape( N, NB, my_pe, src_pe, np, pshape)
-            RETURN
-          END FUNCTION localdim_desc
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-         INTEGER FUNCTION  localindex_shape(ig, n, nb, me, np, pshape)
-
-!  ig  global index of the x dimension of array element
-!  n   dimension of the global array
-!  nb  dimension of the block the global array is split into.
-!  np  number of processors onto which the array is distributed
-!
-!  This function return the index of the element in the local block
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-           INTEGER ig, n, np, pshape, nb, me, q, r
-
-           IF( pshape .EQ. BLOCK_PARTITION_SHAPE ) THEN
-
-             q = INT(n/np)
-             r = MOD(n,np)
-             IF( me < r ) THEN
-               LOCALINDEX_SHAPE = ig - (q+1) * me
-             ELSE
-               LOCALINDEX_SHAPE = ig - (q+1) * r - q * (me - r)
-             END IF
-
-           ELSE IF ( pshape .EQ. BLOCK_CYCLIC_SHAPE ) THEN
-
-             LOCALINDEX_SHAPE = NB*((IG-1)/(NB*NP))+MOD(IG-1,NB)+1
-
-           ELSE IF ( pshape .EQ. CYCLIC_SHAPE ) THEN
- 
-             LOCALINDEX_SHAPE = (ig-1)/np + 1
-
-           ELSE
-             LOCALINDEX_SHAPE = ig
-           END IF
-           RETURN
-         END FUNCTION localindex_shape
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          INTEGER FUNCTION localindex_desc(ig, desc, what )
-
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-            TYPE (descriptor) :: desc
-            CHARACTER(LEN=*) :: what
-            INTEGER ig, n, nb, np, pshape, me
-            IF ( what(1:1) .EQ. 'R' .OR. what(1:1) .EQ. 'r' ) THEN
-              NB = desc%nxblk;      N      = desc%nx; 
-              np = desc%grid%npx;    pshape = desc%xshape
-              me = desc%grid%mex
-            ELSE IF ( what(1:1) .EQ. 'C' .OR. what(1:1) .EQ. 'c' ) THEN
-              NB = desc%nyblk;   N      = desc%ny; 
-              np = desc%grid%npy; pshape = desc%yshape
-              me = desc%grid%mey
-            ELSE IF ( what(1:1) .EQ. 'P' .OR. what(1:1) .EQ. 'p' ) THEN
-              NB = desc%nzblk;    N      = desc%nz; 
-              np = desc%grid%npz;  pshape = desc%zshape
-              me = desc%grid%mez
-            END IF
-            localindex_desc = localindex_shape(ig,n,nb,me,np,pshape)
-            RETURN
-          END FUNCTION localindex_desc
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-         INTEGER FUNCTION ownerof_shape(ig,n,nb,src_pe,np,pshape)
-!
-!  ig      global index of the x dimension of array element
-!  n       dimension of the global array
-!  nb      dimension of the block 
-!  src_pe  index of the processor owning the first element of the array
-!             at the moment meaningfull only for pshape = BLOCK_CYCLIC_SHAPE
-!  np      number of processors
-!
-!  This function return the index of the processor owning the array element
-!  whose global index is "ig"
-!
-!  END manual
-!=----------------------------------------------------------------------------=!
- 
-       IMPLICIT NONE
-       INTEGER ig, n, nb, np, pshape, src_pe, r, q
-       IF( pshape .EQ. BLOCK_PARTITION_SHAPE ) THEN
-         q = INT(n/np);  r = MOD(n,np)
-         IF ( ig <= ((q+1)*r) ) THEN
-           ownerof_shape = INT((ig-1)/(q+1))
-         ELSE
-           ownerof_shape = INT((ig-1-r*(q+1))/q)+r
-         END IF
-       ELSE IF( pshape .EQ. BLOCK_CYCLIC_SHAPE ) THEN
-         ownerof_shape = MOD( src_pe + (ig - 1) / NB, NP )
-       ELSE IF( pshape .EQ. CYCLIC_SHAPE ) THEN
-         ownerof_shape = MOD( ig-1, np )
-       END IF
-       RETURN
-       END FUNCTION ownerof_shape
-
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          INTEGER FUNCTION ownerof_desc(ig, desc, what )
-
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-            TYPE (descriptor) :: desc
-            CHARACTER(LEN=*) :: what
-            INTEGER ig, n, nb, src_pe, np, pshape
-            IF ( what(1:1) .EQ. 'R' .OR. what(1:1) .EQ. 'r' ) THEN
-              NB = desc%nxblk;      N      = desc%nx; 
-              np = desc%grid%npx;    pshape = desc%xshape
-              src_pe = desc%ipexs
-            ELSE IF ( what(1:1) .EQ. 'C' .OR. what(1:1) .EQ. 'c' ) THEN
-              NB = desc%nyblk;   N      = desc%ny; 
-              np = desc%grid%npy; pshape = desc%yshape
-              src_pe = desc%ipeys
-            ELSE IF ( what(1:1) .EQ. 'P' .OR. what(1:1) .EQ. 'p' ) THEN
-              NB = desc%nzblk;    N      = desc%nz; 
-              np = desc%grid%npz;  pshape = desc%zshape
-              src_pe = desc%ipezs
-            END IF
-            ownerof_desc = ownerof_shape(ig, n, nb, src_pe, np, pshape)
-            RETURN
-          END FUNCTION ownerof_desc
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          SUBROUTINE desc_gdims(d, nx, ny, nz )
-
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-            TYPE (descriptor), INTENT(IN) :: d
-            INTEGER, INTENT(OUT) :: nx, ny, nz
-            nx = d%nx
-            ny = d%ny
-            nz = d%nz
-            RETURN
-          END SUBROUTINE desc_gdims
-
-!=----------------------------------------------------------------------------=!
-!  BEGIN manual
-
-          SUBROUTINE desc_ldims(d, nxl, nyl, nzl )
-
-!  END manual
-!=----------------------------------------------------------------------------=!
-
-            TYPE (descriptor), INTENT(IN) :: d
-            INTEGER, INTENT(OUT) :: nxl
-            INTEGER, OPTIONAL, INTENT(OUT) :: nyl, nzl
-            nxl = d%nxl
-            IF( PRESENT( nyl ) ) nyl = d%nyl
-            IF( PRESENT( nzl ) ) nzl = d%nzl
-
-            RETURN
-          END SUBROUTINE desc_ldims
-
-
 
       END MODULE descriptors_module

Modules/electrons_base.f90

@@ -31,7 +31,7 @@
 
       REAL(DP), ALLOCATABLE :: f(:)   ! occupation numbers ( at gamma )
       REAL(DP) :: qbac = 0.0d0        ! background neutralizing charge
-      INTEGER, ALLOCATABLE :: fspin(:) ! spin of each state
+      INTEGER, ALLOCATABLE :: ispin(:) ! spin of each state
 !
 !------------------------------------------------------------------------------!
   CONTAINS
@@ -113,9 +113,9 @@
       END IF
 
       ALLOCATE( f     ( nbspx ) )
-      ALLOCATE( fspin ( nbspx ) )
+      ALLOCATE( ispin ( nbspx ) )
       f     = 0.0d0
-      fspin = 0
+      ispin = 0
 
       iupdwn ( 1 ) = 1
       nel          = 0
@@ -265,7 +265,7 @@
 
       do iss = 1, nspin
          do in = iupdwn(iss), iupdwn(iss) - 1 + nupdwn(iss)
-            fspin(in) = iss
+            ispin(in) = iss
          end do
       end do
 
@@ -402,7 +402,7 @@
 
     SUBROUTINE deallocate_elct()
       IF( ALLOCATED( f ) ) DEALLOCATE( f )
-      IF( ALLOCATED( fspin ) ) DEALLOCATE( fspin )
+      IF( ALLOCATED( ispin ) ) DEALLOCATE( ispin )
       telectrons_base_initval = .FALSE.
       RETURN
     END SUBROUTINE deallocate_elct

Modules/parallel_types.f90

@@ -27,12 +27,13 @@
                               !  0 <= mez < npz-1
         END TYPE
 
-! ...   Valid values for data shape
+        ! ...   Valid values for data distribution
-        INTEGER, PARAMETER :: BLOCK_CYCLIC_SHAPE = 1
+        !
-        INTEGER, PARAMETER :: BLOCK_PARTITION_SHAPE = 2
+        INTEGER, PARAMETER :: BLOCK_CYCLIC_DIST    = 1
-        INTEGER, PARAMETER :: FREE_PATTERN_SHAPE = 3
+        INTEGER, PARAMETER :: BLOCK_PARTITION_DIST = 2
-        INTEGER, PARAMETER :: REPLICATED_DATA_SHAPE = 4
+        INTEGER, PARAMETER :: FREE_PATTERN_DIST    = 3
-        INTEGER, PARAMETER :: CYCLIC_SHAPE = 5
+        INTEGER, PARAMETER :: REPLICATED_DATA_DIST = 4
+        INTEGER, PARAMETER :: CYCLIC_DIST          = 5
 
 !  ----------------------------------------------
 !  BEGIN manual
@@ -40,18 +41,18 @@
 !  Given the Array  |a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11|
 !  and three processors P0, P1, P2 
 !
-!  in the BLOCK_PARTITION_SHAPE scheme, the Array is partitioned 
+!  in the BLOCK_PARTITION_DIST scheme, the Array is partitioned 
 !  as follow
 !       P0            P1            P2
 !  |a1 a2 a3 a4| |a5 a6 a7 a8| |a9 a10 a11|
 !
-!  in the BLOCK_CYCLIC_SHAPE scheme the Array is first partitioned 
+!  in the BLOCK_CYCLIC_DIST scheme the Array is first partitioned 
 !  into blocks (i.e. of size 2)  |a1 a2|a3 a4|a5 a6|a7 a8|a9 a10|a11|
 !  Then the block are distributed cyclically among P0, P1 and P2
 !       P0             P1              P2
 !  |a1 a2|a7 a8|  |a3 a4|a9 a10|  |a5 a6|a11|
 !
-!  in the CYCLIC_SHAPE scheme the Array elements are distributed round robin
+!  in the CYCLIC_DIST scheme the Array elements are distributed round robin
 !  among P0, P1 and P2
 !       P0             P1              P2
 !  |a1 a4 a7 a10|  |a2 a5 a8 a11|  |a3 a6 a9|
@@ -67,7 +68,7 @@
           INTEGER :: nx     ! rows, number of rows in the global array
           INTEGER :: ny     ! columns, number of columns in the global array
           INTEGER :: nz     ! planes, number of planes in the global array
-          INTEGER :: nxblk  ! row_block, if shape = BLOCK_CICLYC_SHAPE,
+          INTEGER :: nxblk  ! row_block, if DIST = BLOCK_CICLYC_DIST,
                             ! this value represent the blocking factor
                             ! used to distribute the rows of the array,
                             ! otherwise this is the size of local block of rows
@@ -90,9 +91,9 @@
                             ! of the array
           INTEGER :: ldz    ! 
 
-          INTEGER :: xshape ! row_shape
+          INTEGER :: xdist ! row_dist
-          INTEGER :: yshape ! column_shape
+          INTEGER :: ydist ! column_dist
-          INTEGER :: zshape ! plane_shape
+          INTEGER :: zdist ! plane_dist
         END TYPE
         
         
@@ -150,8 +151,8 @@
           complex_parallel_vector, complex_parallel_matrix, &
           complex_parallel_tensor, parallel_allocate, parallel_deallocate
 
-        PUBLIC ::  BLOCK_CYCLIC_SHAPE, BLOCK_PARTITION_SHAPE, &
+        PUBLIC ::  BLOCK_CYCLIC_DIST, BLOCK_PARTITION_DIST, &
-          FREE_PATTERN_SHAPE, REPLICATED_DATA_SHAPE, CYCLIC_SHAPE
+          FREE_PATTERN_DIST, REPLICATED_DATA_DIST, CYCLIC_DIST
 
         INTERFACE parallel_allocate
           MODULE PROCEDURE allocate_real_vector, allocate_real_matrix, &

Modules/read_namelists.f90

@@ -1062,8 +1062,8 @@ MODULE read_namelists_module
        IF( calculation == ' ' ) &
           CALL errore( sub_name,' calculation not specified ',1)
        IF( prog == 'CP' ) THEN
-          IF( calculation == 'nscf' .OR. calculation == 'phonon' ) &
+          IF( calculation == 'phonon' ) &
-             CALL errore( sub_name,' calculation '//TRIM(calculation)// &
+             CALL errore( sub_name,' calculation '//calculation// &
                           & ' not implemented ',1)
        END IF
        IF( ndr < 50 ) &
@@ -1573,9 +1573,9 @@ MODULE read_namelists_module
                  cell_dynamics     = 'none'
              END IF
           CASE ('nscf')
-             IF( prog == 'CP' ) &
+             ! IF( prog == 'CP' ) &
-                CALL errore( sub_name,' calculation '//TRIM(calculation)// &
+             !    CALL errore( sub_name,' calculation '//calculation// &
-                             & ' not implemented ',1)
+             !                 & ' not implemented ',1)
              IF( prog == 'CP' ) occupations = 'bogus'
              IF( prog == 'CP' ) electron_dynamics = 'damp'
              IF( prog == 'PW' ) startingpot = 'file'

Modules/xml_io_base.f90

@@ -1001,7 +1001,7 @@ MODULE xml_io_base
       CALL mp_bcast( ierr, ionode_id )
       !
       CALL errore( ' read_rho_xml ', &
-                   'cannot open rho_file file for writing', ierr )
+                   'cannot open ' // rho_file // ' file for reading', ierr )
       !
       IF ( ionode ) THEN
          !

flib/Makefile

@@ -48,7 +48,6 @@ infog1l.o \
 infog2l.o \
 localdim.o \
 localindex.o \
-matmul.o \
 npreroc.o \
 numroc.o \
 ownerof.o \

flib/gridsetup.f90

@@ -10,31 +10,40 @@
 !-----------------------------------------------------------------------
 !
 
-      subroutine GRIDSETUP(NPROC,NPROW,NPCOL)
+SUBROUTINE GRID2D_DIMS( nproc, nprow, npcol )
-
    !
    ! This subroutine factorizes the number of processors (NPROC)
    ! into NPROW and NPCOL,  that are the sizes of the 2D processors mesh.
    !
    !    Written by Carlo Cavazzoni
    !
-
    IMPLICIT NONE
-
+   INTEGER, INTENT(IN)  :: nproc
-      integer nproc,nprow,npcol
+   INTEGER, INTENT(OUT) :: nprow, npcol
    integer sqrtnp,i
-
+   sqrtnp = int( sqrt( dble(nproc) ) + 1 )
-      if(nproc.lt.2) then
-        npcol = 1
-        nprow = 1
-      else
-        sqrtnp = int( sqrt( DBLE(nproc) ) + 1 )
    do i=1,sqrtnp
       if(mod(nproc,i).eq.0) nprow = i
    end do
    npcol = nproc/nprow
-      endif
+   RETURN
+END SUBROUTINE
 
-      return
+SUBROUTINE GRID2D_COORDS( rank, nprow, npcol, row, col )
-      end subroutine gridsetup
+   IMPLICIT NONE
+   INTEGER, INTENT(IN)  ::  rank          ! process index starting from 0
+   INTEGER, INTENT(IN)  ::  nprow, npcol  ! dimensions of the processor grid
+   INTEGER, INTENT(OUT) ::  row, col
+   row = MOD( rank, nprow )
+   col = rank / nprow
+   RETURN
+END SUBROUTINE
 
+SUBROUTINE GRID2D_RANK( nprow, npcol, row, col, rank )
+   IMPLICIT NONE
+   INTEGER, INTENT(OUT) ::  rank         ! process index starting from 0
+   INTEGER, INTENT(IN)  ::  nprow, npcol  ! dimensions of the processor grid
+   INTEGER, INTENT(IN)  ::  row, col
+   rank = row + col * nprow
+   RETURN
+END SUBROUTINE

flib/localindex.f90

@@ -122,3 +122,83 @@
 !     End of INDXG2L
 !
       END FUNCTION lind_block_cyclic
+
+
+!=----------------------------------------------------------------------------=!
+
+
+      INTEGER FUNCTION gind_cyclic( lind, n, np, me )
+
+!  This function computes the global index of a distributed array entry
+!  pointed to by the local index lind of the process indicated by me.
+!  lind      local index of the distributed matrix entry.
+!  N         is the size of the global array.
+!  me        The coordinate of the process whose local array row or
+!            column is to be determined.
+!  np        The total number processes over which the distributed
+!            matrix is distributed.
+!
+
+            INTEGER, INTENT(IN) :: lind, n, me, np
+            INTEGER r, q
+
+            gind_cyclic = (lind-1) * np + me + 1
+
+            RETURN
+      END FUNCTION gind_cyclic
+
+
+!=----------------------------------------------------------------------------=!
+
+
+      INTEGER FUNCTION gind_block( lind, n, np, me )
+
+!  This function computes the global index of a distributed array entry
+!  pointed to by the local index lind of the process indicated by me.
+!  lind      local index of the distributed matrix entry.
+!  N         is the size of the global array.
+!  me        The coordinate of the process whose local array row or
+!            column is to be determined.
+!  np        The total number processes over which the distributed
+!            matrix is distributed.
+
+
+            INTEGER, INTENT(IN) :: lind, n, me, np
+            INTEGER r, q
+
+              q = INT(n/np)
+              r = MOD(n,np)
+              IF( me < r ) THEN
+                gind_block = (Q+1)*me + lind
+              ELSE
+                gind_block = Q*me + R + lind
+              END IF
+
+         RETURN
+      END FUNCTION gind_block
+
+!=----------------------------------------------------------------------------=!
+
+      INTEGER FUNCTION gind_block_cyclic( lind, n, nb, np, me )
+
+!  This function computes the global index of a distributed array entry
+!  pointed to by the local index lind of the process indicated by me.
+!  lind      local index of the distributed matrix entry.
+!  N         is the size of the global array.
+!  NB        size of the blocks the distributed matrix is split into.
+!  me        The coordinate of the process whose local array row or
+!            column is to be determined.
+!  np        The total number processes over which the distributed
+!            matrix is distributed.
+
+
+            INTEGER, INTENT(IN) :: lind, n, nb, me, np
+            INTEGER r, q, isrc
+
+            isrc = 0
+            gind_block_cyclic = np*NB*((lind-1)/NB) + &
+                MOD(lind-1,NB) + MOD(np+me-isrc, np)*NB + 1
+
+         RETURN
+      END FUNCTION gind_block_cyclic
+

pwtools/matdyn.f90

@@ -1614,7 +1614,8 @@ SUBROUTINE a2Fdos &
   !
   INTEGER, INTENT(in) :: nat, nq, nr1, nr2, nr3, ibrav, ndos, ntetra, &
        tetra(4, ntetra)
-  LOGICAL, INTENT(in) :: dos, asr
+  LOGICAL, INTENT(in) :: dos
+  CHARACTER(LEN=*), INTENT(IN) :: asr
   REAL(DP), INTENT(in) :: freq(3*nat,nq), q(3,nq), at(3,3), bg(3,3), &
        tau(3,nat), alat, Emin, DeltaE
   !