Reintroduced the changes to the paw pseudopotential generetion, not merged

in the commit of two days ago.


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

Modules/parameters.f90

@@ -21,7 +21,8 @@ MODULE parameters
        npk    = 40000,  &! max number of k-points               
        lmaxx  = 3,      &! max non local angular momentum (l=0 to lmaxx)      
        cp_lmax  = lmaxx+1,&! max number of channels !TMP FOR PAW: REMOVE
-       nchix  = 6        ! max number of atomic wavefunctions per atom
+       nchix  = 6,      & ! max number of atomic wavefunctions per atom
+       nwfsx  = 14       ! max number of beta functions per atom
   !
   INTEGER, PARAMETER :: &
        nbrx = 14,          &! max number of beta functions

Modules/pseudo_types.f90

@@ -17,8 +17,9 @@
         USE kinds, ONLY: DP
         USE parameters, ONLY: cp_lmax, lmaxx
         use radial_grids, ONLY: ndmx, radial_grid_type
-        !  USE ld1_parameters, ONLY: ndm, nwfsx
+!        USE ld1_parameters, ONLY: nwfsx
-        USE parameters, ONLY: nwfsx=>nchix
+        USE parameters, ONLY: nwfsx
+!        USE ld1_parameters, ONLY: nwfsx
 
         IMPLICIT NONE
         SAVE

Modules/read_paw.f90

@@ -65,6 +65,12 @@
        WRITE(un,'(i8)')     (pawset_%ikk(ns), ns=1,pawset_%nwfc)
         write(un,'(a)') "oc:"
        WRITE(un,'(e20.10)') (pawset_%oc(ns), ns=1,pawset_%nwfc)
+        write(un,'(a)') "els:"
+       WRITE(un,'(a2)') (pawset_%els(ns), ns=1,pawset_%nwfc)
+        write(un,'(a)') "jj:"
+       WRITE(un,'(e20.10)') (pawset_%jj(ns), ns=1,pawset_%nwfc)
+        write(un,'(a)') "rcutus:"
+       WRITE(un,'(e20.10)') (pawset_%rcutus(ns), ns=1,pawset_%nwfc)
         write(un,'(a)') "enl:"
        WRITE(un,'(e20.10)') (pawset_%enl(ns), ns=1,pawset_%nwfc)
         write(un,'(a)') "aewfc:"
@@ -145,6 +151,12 @@
        READ (un,'(i8)')     (pawset_%ikk(ns), ns=1,pawset_%nwfc)
         read(un, '(a)') dummy
        READ (un,'(e20.10)') (pawset_%oc(ns), ns=1,pawset_%nwfc)
+        read(un,'(a)') dummy
+       READ(un,'(a2)') (pawset_%els(ns), ns=1,pawset_%nwfc)
+        read(un,'(a)') dummy
+       READ(un,'(e20.10)') (pawset_%jj(ns), ns=1,pawset_%nwfc)
+        read(un,'(a)') dummy
+       READ(un,'(e20.10)') (pawset_%rcutus(ns), ns=1,pawset_%nwfc)
         read(un, '(a)') dummy
        READ (un,'(e20.10)') (pawset_%enl(ns), ns=1,pawset_%nwfc)
         read(un, '(a)') dummy
@@ -247,8 +259,11 @@ SUBROUTINE allocate_pseudo_paw( paw, size_mesh, size_nwfc, size_lmax )
   INTEGER, INTENT(IN) :: size_mesh, size_nwfc, size_lmax
 !   WRITE(0,"(a,3i5)") "Allocating PAW setup: ",size_mesh, size_nwfc, size_lmax
   ALLOCATE ( paw%l(size_nwfc) )
+  ALLOCATE ( paw%jj(size_nwfc) )
   ALLOCATE ( paw%ikk(size_nwfc) )
   ALLOCATE ( paw%oc(size_nwfc) )
+  ALLOCATE ( paw%rcutus(size_nwfc) )
+  ALLOCATE ( paw%els(size_nwfc) )
   ALLOCATE ( paw%enl(size_nwfc) )
   ALLOCATE ( paw%aewfc(size_mesh,size_nwfc) )
   ALLOCATE ( paw%pswfc(size_mesh,size_nwfc) )
@@ -269,8 +284,11 @@ SUBROUTINE deallocate_pseudo_paw( paw )
 
   TYPE( paw_t ), INTENT(INOUT) :: paw
   IF( ASSOCIATED( paw%l ) ) DEALLOCATE( paw%l )
+  IF( ASSOCIATED( paw%jj ) ) DEALLOCATE( paw%jj )
   IF( ASSOCIATED( paw%ikk ) ) DEALLOCATE( paw%ikk )
   IF( ASSOCIATED( paw%oc ) ) DEALLOCATE( paw%oc )
+  IF( ASSOCIATED( paw%els ) ) DEALLOCATE( paw%els )
+  IF( ASSOCIATED( paw%rcutus ) ) DEALLOCATE( paw%rcutus )
   IF( ASSOCIATED( paw%enl ) ) DEALLOCATE( paw%enl )
   IF( ASSOCIATED( paw%aewfc ) ) DEALLOCATE( paw%aewfc )
   IF( ASSOCIATED( paw%pswfc ) ) DEALLOCATE( paw%pswfc )

atomic/atomic_paw.f90

@@ -69,7 +69,7 @@ CONTAINS
   ! Compute also the total energy
   ! 
   SUBROUTINE new_paw_hamiltonian (veffps_, ddd_, etot_, &
-       pawset_, nwfc_, l_, nspin_, spin_, oc_, pswfc_, eig_, dddion_)
+       pawset_, nwfc_, l_, nspin_, spin_, oc_, pswfc_, eig_, paw_energy,dddion_)
     IMPLICIT NONE
     REAL(dp), INTENT(OUT) :: veffps_(ndmx,2)
     REAL(dp), INTENT(OUT) :: ddd_(nwfsx,nwfsx,2)
@@ -83,6 +83,7 @@ CONTAINS
     REAL(dp), INTENT(IN)  :: pswfc_(ndmx,nwfsx)
     REAL(dp), INTENT(IN)  :: eig_(nwfsx)
     REAL(dp), OPTIONAL :: dddion_(nwfsx,nwfsx,2)
+    REAL(dp), INTENT(OUT), OPTIONAL :: paw_energy(5,3) 
     !
     REAL(dp) :: &                                        ! one center:
          eps,             e1,             e1ps,             & ! energies;
@@ -91,11 +92,11 @@ CONTAINS
          projsum(nwfsx,nwfsx,2), eigsum !  sum of projections, sum of eigenval.
     !
     INTEGER :: ns, ns1, is
-    REAL(dp) :: aux(ndmx), energies(4)
+    REAL(dp) :: aux(ndmx), energy(5,3)
     !
     ! Compute the valence charges
     CALL compute_charges(projsum, chargeps, charge1, charge1ps, &
-       pawset_, nwfc_, l_, nspin_, spin_, oc_, pswfc_ )
+       pawset_, nwfc_, l_, nspin_, spin_, oc_, pswfc_, 1 )
 !     write(766,"(4f12.6)") (pawset_%grid%r(ns), chargeps(ns,1), charge1(ns,1), charge1ps(ns,1), ns=1,pawset_%grid%mesh)
 !     write(767,"(4f12.6)") (pawset_%grid%r(ns), chargeps(ns,2), charge1(ns,2), charge1ps(ns,2), ns=1,pawset_%grid%mesh)
     !
@@ -105,13 +106,16 @@ CONTAINS
     ! where n_v can be n~+n^ or n1 or n1~+n^ and n_c can be nc or nc~
     ! n~+n^ , nc~
     CALL compute_onecenter_energy ( eps,  veffps_, &
-       pawset_, chargeps,  pawset_%nlcc, pawset_%psccharge, nspin_)
+       pawset_, chargeps,  pawset_%nlcc, pawset_%psccharge, nspin_,&
+             pawset_%grid%mesh, pawset_%psloc, energy(1,1))
     ! n1 , nc
     CALL compute_onecenter_energy ( e1,   veff1, &
-       pawset_, charge1,  .TRUE.,        pawset_%aeccharge, nspin_)
+       pawset_, charge1,  .TRUE.,        pawset_%aeccharge, nspin_,&
+            pawset_%irc, pawset_%aeloc, energy(1,2))
     ! n1~+n^ , nc~
     CALL compute_onecenter_energy ( e1ps, veff1ps, &
-       pawset_, charge1ps, pawset_%nlcc, pawset_%psccharge, nspin_)
+       pawset_, charge1ps, pawset_%nlcc, pawset_%psccharge, nspin_,&
+            pawset_%irc, pawset_%psloc, energy(1,3))
     ! Add the local part
     DO is=1,nspin_
        veffps_(1:pawset_%grid%mesh,is) = veffps_(1:pawset_%grid%mesh,is) +  &
@@ -134,6 +138,8 @@ CONTAINS
        IF (oc_(ns)>ZERO) eigsum = eigsum + oc_(ns)*eig_(ns)
     END DO
     etot_ = eigsum + eps + e1 - e1ps
+
+    if (PRESENT(paw_energy)) paw_energy=energy
     !
   END SUBROUTINE new_paw_hamiltonian
   !
@@ -143,7 +149,8 @@ CONTAINS
   !
   SUBROUTINE us2paw (pawset_,                                     &
        zval, grid, rmatch_augfun, ikk, &
-       nbeta, lls, ocs, enls, psipaw, phis, betas, qvan, kindiff, &
+       nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, phis, betas, &
+       qvan, kindiff, &
        nlcc, aerhoc, psrhoc, aevtot, psvtot, which_paw_augfun)
     USE funct, ONLY : dft_name, get_iexch, get_icorr, get_igcx, get_igcc
     USE ld1inc, ONLY : zed
@@ -159,6 +166,9 @@ CONTAINS
     INTEGER,  INTENT(IN)  :: nbeta
     INTEGER,  INTENT(IN)  :: lls(nwfsx)
     REAL(dp), INTENT(IN)  :: ocs(nwfsx)
+    CHARACTER(LEN=2), INTENT(IN)  :: els(nwfsx)
+    REAL(dp), INTENT(IN)  :: jjs(nwfsx)
+    REAL(dp), INTENT(IN)  :: rcutus(nwfsx)
     REAL(dp), INTENT(IN)  :: enls(nwfsx)
     REAL(dp), INTENT(IN)  :: psipaw(ndmx,nwfsx)
     REAL(dp), INTENT(IN)  :: phis(ndmx,nwfsx)
@@ -184,7 +194,7 @@ CONTAINS
     !
     ! variables for aug. functions generation
     ! 
-    REAL(DP) :: qc(2), xc(2), b1(2), b2(2)
+    REAL(DP) :: qc(2), xc(2), b1(2), b2(2), energy(5,3)
     REAL(DP), ALLOCATABLE :: j1(:,:)
     INTEGER  :: nc, iok          ! index Bessel function, ...
     INTEGER :: l1,l2,l3, lll, ircm, ir, ir0
@@ -218,6 +228,9 @@ CONTAINS
     pawset_%l(1:nbeta) = lls(1:nbeta)
     pawset_%lmax = MAXVAL(pawset_%l(1:pawset_%nwfc))
     pawset_%oc(1:nbeta) = ocs(1:nbeta)
+    pawset_%jj(1:nbeta) = jjs(1:nbeta)
+    pawset_%rcutus(1:nbeta) = rcutus(1:nbeta)
+    pawset_%els(1:nbeta)= els(1:nbeta)
     pawset_%enl(1:nbeta)= enls(1:nbeta)
     pawset_%aewfc(1:mesh,1:nbeta) = psipaw(1:mesh,1:nbeta)
     pawset_%pswfc(1:mesh,1:nbeta) = phis  (1:mesh,1:nbeta)
@@ -378,10 +391,12 @@ CONTAINS
        pawset_, nbeta, lls, nspin, spin, ocs, phis )
     pawset_%pscharge(1:mesh)=pscharge(1:mesh,1)
     CALL compute_onecenter_energy ( raux,  aux2, &
-       pawset_, pscharge, pawset_%nlcc, pawset_%psccharge, nspin )
+       pawset_, pscharge, pawset_%nlcc, pawset_%psccharge, nspin, &
+       pawset_%grid%mesh )
     pawset_%psloc(1:mesh)=psvtot(1:mesh)-aux2(1:mesh,1)
     CALL compute_onecenter_energy ( raux,  aux2, &
-       pawset_, aecharge, .TRUE.,       pawset_%aeccharge, nspin )
+       pawset_, aecharge, .TRUE.,       pawset_%aeccharge, nspin, &
+       pawset_%grid%mesh)
     pawset_%aeloc(1:mesh)=aevtot(1:mesh)-aux2(1:mesh,1)
     !WRITE(4444,'(5e20.10)')(r(n),aevtot(n),psvtot(n),pawset_%aeloc(n),pawset_%psloc(n),n=1,mesh)
     !
@@ -394,7 +409,7 @@ CONTAINS
     !
     ! Generate the paw hamiltonian for test (should be equal to the US one)
     CALL new_paw_hamiltonian (vps, ddd, etot, &
-       pawset_, pawset_%nwfc, pawset_%l, nspin, spin, pawset_%oc, pawset_%pswfc, pawset_%enl, dddion)
+       pawset_, pawset_%nwfc, pawset_%l, nspin, spin, pawset_%oc, pawset_%pswfc, pawset_%enl, energy, dddion)
     pawset_%dion(1:nbeta,1:nbeta)=dddion(1:nbeta,1:nbeta,1)
     WRITE(6,'(/5x,A,f12.6,A)') 'Estimated PAW energy =',etot,' Ryd'
     WRITE(6,'(/5x,A)') 'The PAW screened D coefficients'
@@ -413,8 +428,8 @@ CONTAINS
   !
   ! ...
   !
-  SUBROUTINE paw2us (pawset_,zval,grid,nbeta,lls,ikk,betas,qq,qvan,&
+  SUBROUTINE paw2us (pawset_,zval,grid,nbeta,lls,jjs,ikk,betas,qq,qvan,&
-                     vpsloc, bmat, rhos, pseudotype)
+                     vpsloc, bmat, rhos, els, rcutus, pseudotype)
     USE funct, ONLY : set_dft_from_name
     use radial_grids, only: radial_grid_type
     IMPLICIT NONE
@@ -425,6 +440,9 @@ CONTAINS
     INTEGER,       INTENT(OUT) :: lls(nwfsx)
     INTEGER,       INTENT(OUT) :: ikk(nwfsx)
     INTEGER,       INTENT(OUT) :: pseudotype
+    CHARACTER(LEN=2), INTENT(OUT) :: els(nwfsx)
+    REAL(dp), INTENT(OUT) :: jjs(nwfsx)
+    REAL(dp), INTENT(OUT) :: rcutus(nwfsx)
     REAL(dp), INTENT(OUT) :: betas(ndmx,nwfsx)
     REAL(dp), INTENT(OUT) :: qq(nwfsx,nwfsx)
     REAL(dp), INTENT(OUT) :: qvan(ndmx,nwfsx,nwfsx)
@@ -454,6 +472,9 @@ CONTAINS
     !
     nbeta=pawset_%nwfc
     lls(1:nbeta)=pawset_%l(1:nbeta)
+    jjs(1:nbeta)=pawset_%jj(1:nbeta)
+    els(1:nbeta)=pawset_%els(1:nbeta)
+    rcutus(1:nbeta)=pawset_%rcutus(1:nbeta)
     ikk(1:nbeta)=pawset_%ikk(1:nbeta)
     !
     DO ns=1,nbeta
@@ -536,7 +557,7 @@ CONTAINS
   !============================================================================
   !
   SUBROUTINE compute_charges (projsum_, chargeps_, charge1_, charge1ps_, &
-       pawset_, nwfc_, l_, nspin_, spin_, oc_, pswfc_ , unit_)
+       pawset_, nwfc_, l_, nspin_, spin_, oc_, pswfc_ , iflag, unit_)
     IMPLICIT NONE
     REAL(dp), INTENT(OUT) :: projsum_(nwfsx,nwfsx,2)
     REAL(dp), INTENT(OUT) :: chargeps_(ndmx,2)
@@ -549,12 +570,15 @@ CONTAINS
     INTEGER,       INTENT(IN)  :: spin_(nwfsx)
     REAL(dp), INTENT(IN)  :: oc_(nwfsx)
     REAL(dp), INTENT(IN)  :: pswfc_(ndmx,nwfsx)
-    INTEGER, OPTIONAL :: unit_
+    INTEGER, OPTIONAL :: unit_, iflag
-    REAL(dp) :: augcharge(ndmx,2)
+    REAL(dp) :: augcharge(ndmx,2), chargetot
-    INTEGER :: i
+    INTEGER :: i, n, ns, ns1, iflag0
+
+    iflag0=0
+    if (present(iflag)) iflag0=iflag
     CALL compute_sumwfc2(chargeps_,pawset_,nwfc_,pswfc_,oc_,spin_)
     CALL compute_projsum(projsum_,pawset_,nwfc_,l_,spin_,pswfc_,oc_)
-    !WRITE (6200,'(20e20.10)') ((projsum(ns,ns1),ns=1,ns1),ns1=1,pawset_%nwfc)
+!    WRITE (6200,'(20e20.10)') ((projsum_(ns,ns1,1),ns=1,ns1),ns1=1,pawset_%nwfc)
     CALL compute_onecenter_charge(charge1ps_,pawset_,projsum_,nspin_,"PS")
     CALL compute_onecenter_charge(charge1_  ,pawset_,projsum_,nspin_,"AE")
     ! add augmentation charges
@@ -571,6 +595,26 @@ CONTAINS
                                         + augcharge(1:pawset_%grid%mesh,1:nspin_)
     charge1ps_(1:pawset_%grid%mesh,1:nspin_) = charge1ps_(1:pawset_%grid%mesh,1:nspin_) &
                                         + augcharge(1:pawset_%grid%mesh,1:nspin_)
+!
+!  If there are unbounded scattering states in the pseudopotential generation,
+!  n1 and n~1 separately diverge. This makes the hartree and exchange and
+!  correlation energies separately diverging. The cancellation becomes
+!  very difficult numerically. Outside the sphere however the two charges 
+!  should be equal and opposite and we set them to zero. 
+!  Is this the right solution? 
+!
+    do n=pawset_%irc+1,pawset_%grid%mesh
+       chargetot=chargeps_(n,1)
+       if (nspin_==2) chargetot=chargetot+chargeps_(n,2)
+       if (chargetot<1.d-11.or.iflag0==1) then
+          charge1_(n,1:nspin_)=0.0_DP
+          charge1ps_(n,1:nspin_)=0.0_DP
+       endif
+    enddo
+!    do n=1,pawset_%grid%mesh
+!       write(6,'(4f20.15)') pawset_%grid%r(n), chargeps_(n,1), charge1_(n,1),
+!                                               charge1ps_(n,1)
+
   END SUBROUTINE compute_charges
   !
   !============================================================================
@@ -581,24 +625,27 @@ CONTAINS
   ! where n_v can be n~+n^ or n1 or n1~+n^ and n_c can be nc or n~c
   !
   SUBROUTINE compute_onecenter_energy ( totenergy_, veff_, &
-       pawset_, vcharge_, nlcc_, ccharge_, nspin_, energies_ , unit_)
+       pawset_, vcharge_, nlcc_, ccharge_, nspin_, iint, vloc, energies_ , unit_)
     USE funct, ONLY: dft_is_gradient !igcx, igcc
     USE radial_grids, ONLY: hartree
     IMPLICIT NONE
-    REAL(dp), INTENT(OUT) :: totenergy_       ! H+XC+DC
+    REAL(dp), INTENT(OUT) :: totenergy_            ! H+XC+DC
-    REAL(dp), INTENT(OUT) :: veff_(ndmx,2)     ! effective potential
+    REAL(dp), INTENT(OUT) :: veff_(ndmx,2)         ! effective potential
     TYPE(paw_t),   INTENT(IN)  :: pawset_          ! the PAW dataset
-    REAL(dp), INTENT(IN)  :: vcharge_(ndmx,2)  ! valence charge
+    REAL(dp), INTENT(IN)  :: vcharge_(ndmx,2)      ! valence charge
     LOGICAL,       INTENT(IN)  :: nlcc_            ! non-linear core correction
-    REAL(dp), INTENT(IN)  :: ccharge_(ndmx)    ! core charge
+    REAL(dp), INTENT(IN)  :: ccharge_(ndmx)        ! core charge
     INTEGER,       INTENT(IN)  :: nspin_           ! 1 for LDA, 2 for LSDA
-    REAL(dp), INTENT(OUT), OPTIONAL :: energies_(4) ! Etot,H,XC,DC terms
+    INTEGER,       INTENT(IN)  :: iint             ! integrals up to iint
+    REAL(dp), INTENT(IN), OPTIONAL :: vloc(ndmx)   ! 
+    REAL(dp), INTENT(OUT), OPTIONAL :: energies_(5)! Etot,H,XC,DC terms
     INTEGER, OPTIONAL :: unit_
     !
     REAL(dp), PARAMETER :: rho_eq_0(ndmx) = ZERO ! ccharge=0 when nlcc=.f.
     !
     REAL(dp) :: &
          eh, exc, edc, & ! hartree, xc and double counting energies
+         eloc,         & ! local energy
          rhovtot(ndmx), & ! total valence charge
          aux(ndmx),     & ! auxiliary to compute integrals
          vh(ndmx),      & ! hartree potential
@@ -642,7 +689,7 @@ CONTAINS
 #endif
     vh(1:pawset_%grid%mesh) = e2 * vh(1:pawset_%grid%mesh)
     aux(1:pawset_%grid%mesh) = vh(1:pawset_%grid%mesh) * rhovtot(1:pawset_%grid%mesh)
-    eh = HALF * int_0_inf_dr(aux,pawset_%grid,pawset_%grid%mesh,2)
+    eh = HALF * int_0_inf_dr(aux,pawset_%grid,iint,2)
     !
     ! Exhange-Correlation
     rh=(/ZERO,ZERO/)
@@ -663,25 +710,37 @@ CONTAINS
     END DO
     IF (dft_is_gradient()) THEN
        IF (nlcc_) THEN
-          CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,pawset_%grid%r2,vcharge_,ccharge_,vgc,egc,0)
+          CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,&
+                     pawset_%grid%r2,vcharge_,ccharge_,vgc,egc,1)
        ELSE
-          CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,pawset_%grid%r2,vcharge_,rho_eq_0,vgc,egc,0)
+          CALL vxcgc(ndmx,pawset_%grid%mesh,nspin_,pawset_%grid%r,&
+                     pawset_%grid%r2,vcharge_,rho_eq_0,vgc,egc,1)
        END IF
        vxc(1:pawset_%grid%mesh,1:nspin_) = vxc(1:pawset_%grid%mesh,1:nspin_) + &
                                       vgc(1:pawset_%grid%mesh,1:nspin_)
        aux(1:pawset_%grid%mesh) = aux(1:pawset_%grid%mesh) + &
            egc(1:pawset_%grid%mesh) * pawset_%grid%r2(1:pawset_%grid%mesh) * FPI
     END IF
-    exc = int_0_inf_dr(aux,pawset_%grid,pawset_%grid%mesh,2)
+    exc = int_0_inf_dr(aux,pawset_%grid,iint,2)
     !
     ! Double counting
     edc=ZERO
     DO is=1,nspin_
        veff_(1:pawset_%grid%mesh,is)=vxc(1:pawset_%grid%mesh,is)+vh(1:pawset_%grid%mesh)
        aux(1:pawset_%grid%mesh)=veff_(1:pawset_%grid%mesh,is)*vcharge_(1:pawset_%grid%mesh,is)
-       edc=edc+int_0_inf_dr(aux,pawset_%grid,pawset_%grid%mesh,2)
+       edc=edc+int_0_inf_dr(aux,pawset_%grid,iint,2)
     END DO
     !
+    eloc=ZERO
+    !
+    IF (present(vloc)) THEN
+       DO is=1,nspin_
+          aux(1:pawset_%grid%mesh)=vloc(1:pawset_%grid%mesh)        &
+                               *vcharge_(1:pawset_%grid%mesh,is)
+          eloc=eloc+int_0_inf_dr(aux,pawset_%grid,iint,2)
+       ENDDO
+    ENDIF
+    !
     ! Total
     totenergy_ = eh + exc - edc
     !
@@ -690,6 +749,7 @@ CONTAINS
        energies_(2)=eh
        energies_(3)=exc
        energies_(4)=edc
+       energies_(5)=eloc
     END IF
     !
   END SUBROUTINE compute_onecenter_energy
@@ -861,7 +921,7 @@ CONTAINS
                 aux(nr)=pawset_%proj(nr,ns)*pswfc_(nr,nf)
              END DO
              nst=(l_(nf)+1)*2
-             proj_dot_wfc(ns,nf)=int_0_inf_dr(aux,pawset_%grid,pawset_%irc,nst)
+             proj_dot_wfc(ns,nf)=int_0_inf_dr(aux,pawset_%grid,pawset_%ikk(ns),nst)
           ELSE
              proj_dot_wfc(ns,nf)=ZERO
           END IF

atomic/elsdps_paw.f90

@@ -12,8 +12,6 @@
 !
 !   total paw energy in the local-spin-density scheme
 !
-call infomsg('elsdps_paw', 'energy report for PAW coming soon...')
-#ifdef __DONT_DEFINE_ME
 use kinds, only: DP
 use constants, only: fpi
 use radial_grids, only : ndmx
@@ -93,6 +91,5 @@ do ns=1,nwfts
       ekin=ekin+octs(ns)*enlts(ns)
    endif
 end do
-#endif
 return
 end subroutine elsdps_paw

atomic/gener_pseudo.f90

@@ -423,8 +423,9 @@ subroutine gener_pseudo
      ! create the 'pawsetup' object containing the atomic setup for PAW
      call us2paw ( pawsetup,                                         &
           zval, grid, paw_rmatch_augfun, ikk,  &
-          nbeta, lls, ocs, enls, psipaw, phis, betas, qvan, kindiff, &
+          nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, phis, betas, &
-          nlcc, aeccharge, psccharge, vpotpaw, vpsloc, which_paw_augfun)
+          qvan, kindiff, nlcc, aeccharge, psccharge, vpotpaw, vpsloc, &
+          which_paw_augfun)
      !
   endif
   !
@@ -432,8 +433,8 @@ subroutine gener_pseudo
   !
   if (lpaw) then
      ! reread augmentation functions and descreened potentials from PAW
-     call paw2us ( pawsetup, zval, grid, nbeta, lls, ikk, betas, &
+     call paw2us ( pawsetup, zval, grid, nbeta, lls, jjs, ikk, betas, &
-                   qq, qvan, vpsloc, bmat, rhos, pseudotype )
+                   qq, qvan, vpsloc, bmat, rhos, els, rcutus, pseudotype )
   else
      call descreening
   end if

atomic/ld1_readin.f90

@@ -533,8 +533,8 @@ subroutine ld1_readin
         !call paw_io(pawsetup,111,"INP")
         call paw_io(pawsetup,111,"INP",ndmx,nwfsx,lmaxx)
         close(111)
-        call paw2us ( pawsetup, zval, grid, nbeta, lls, ikk, betas, &
+        call paw2us ( pawsetup, zval, grid, nbeta, lls, jjs, ikk, betas, &
-                      qq, qvan,vpsloc, bmat, rhos, pseudotype )
+                      qq, qvan,vpsloc, bmat, rhos, els, rcutus, pseudotype )
         call check_mesh(grid)
         !
      else if ( matches('.rrkj3', file_pseudo) .or. &

atomic/ld1_setup.f90

@@ -23,7 +23,7 @@ subroutine ld1_setup
                 nwf,                       jj,   el,   isw,   oc, nstoae, &
                 nwfs,          lls,        jjs,  els,  isws,  ocs, &
                 nwfts,  nnts,  llts, nnts, jjts, elts, iswts, octs, nstoaets, &
-                nwftsc, nntsc, lltsc, jjtsc, eltsc, iswtsc, octsc, nstoaec
+                nwftsc, nntsc, lltsc, jjtsc, eltsc, iswtsc, octsc, nstoaec, lpaw
   use funct, only : get_iexch, dft_is_meta, start_exx !, set_dft_from_name
   implicit none
 
@@ -138,6 +138,9 @@ subroutine ld1_setup
            &            'the wavefunction chosen for local potential',nsloc)
            end if
         end if
+        if (lpaw .and. ocs(nsloc)>0.0_DP) &
+            call errore('ld1_setup','Paw generation with electrons' // &
+            &    'in the local channel is not available',1)
      else
         nsloc=-1   
         nbeta=nwfs

atomic/ld1inc.f90

@@ -248,8 +248,10 @@ module ld1inc
        psipaw(ndmx,nwfsx),& ! the all-electron wavefunctions for any beta
        aeccharge(ndmx),   & ! true, not smoothened, AE core charge for PAW
        psccharge(ndmx),   & ! smoothened core charge for PAW
-       rCutNC2paw(nwfsx)    ! a cut-off radius for NC wavefunctions to be used
+       rCutNC2paw(nwfsx), & ! a cut-off radius for NC wavefunctions to be used
                             ! instead of AE ones in the construction of PAW
+       paw_energy(5,3)
+
    character(len=20) ::&
        which_paw_augfun     ! choose shape of paw aug.fun. (GAUSS, BESSEL..)
   !

atomic/parameters.f90

@@ -8,6 +8,6 @@
 module ld1_parameters
    integer, parameter :: &
            ncmax1=10,    & ! the maximum configuration number
-           nwfsx=10,     & ! the maximum number of pseudo wavefunctions
+           nwfsx=14,     & ! the maximum number of pseudo wavefunctions
-           nwfx=35         ! the maximum number of wavefunctions
+           nwfx=38         ! the maximum number of wavefunctions
 end module ld1_parameters

atomic/run_pseudo.f90

@@ -21,7 +21,8 @@ subroutine run_pseudo
                      nstoaets, grid, nspin, iter, rhos, rhoc, &
                      nwfts, enlts, llts, jjts, iswts, octs, phits, &
                      vxt, enne, vh, vpsloc, file_potscf, beta, tr2,  &
-                     eps0, file_recon, deld, vpstot, nbeta, ddd, etots
+                     eps0, file_recon, deld, vpstot, nbeta, ddd, etots, &
+                     paw_energy
   use atomic_paw, only : new_paw_hamiltonian
   implicit none
 
@@ -44,7 +45,6 @@ subroutine run_pseudo
   real(DP) :: &
        nvalts,                  & ! number of valence electrons for this conf.
        dddnew(nwfsx,nwfsx,2),   & ! the new D coefficients
-       ocstart(nwfsx),          & ! guess for the occupations
        vd(2*(ndmx+nwfsx+nwfsx)), & ! Vloc and D in one array for mixing
        vdnew(2*(ndmx+nwfsx+nwfsx)) ! the new vd array
   integer :: &
@@ -63,24 +63,12 @@ subroutine run_pseudo
   !
   !    compute an initial estimate of the potential
   !
+  call guess_initial_wfc()
   if (.not.lpaw) then
-     call guess_initial_wfc()
      call start_potps ( )
   else
-     ! Set starting occupations by rescaling those of the generating configuration
+     CALL new_paw_hamiltonian (vpstot, ddd, etots,pawsetup, nwfts, &
-     nvalts=0._dp
+                               llts, nspin, iswts, octs, phits, enlts)
-     do ns=1,nbeta
-        if (octs(ns)>0._dp) nvalts=nvalts+octs(ns)
-     end do
-     ocstart(1:pawsetup%nwfc) = pawsetup%oc(1:pawsetup%nwfc) * nvalts &
-                               / SUM(pawsetup%oc(1:pawsetup%nwfc))
-     iswstart=1
-     ! Generate the corresponding local and nonlocal potentials
-     CALL new_paw_hamiltonian (vpstot, ddd, etots, &
-          pawsetup, pawsetup%nwfc, pawsetup%l, 1,iswstart,ocstart, &
-          pawsetup%pswfc, pawsetup%enl)
-     vpstot(1:grid%mesh,2)=vpstot(1:grid%mesh,1)
-     ddd(1:nbeta,1:nbeta,2)=ddd(1:nbeta,1:nbeta,1)
      do is=1,nspin
         vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)-pawsetup%psloc(1:grid%mesh)
      enddo
@@ -174,9 +162,9 @@ subroutine run_pseudo
   if (.not.lpaw) then
      call elsdps ( )
   else
-      call new_paw_hamiltonian (vnew, dddnew, etots, &
+     call new_paw_hamiltonian (vnew, dddnew, etots, pawsetup, nwfts, &
-           pawsetup, nwfts, llts, nspin, iswts, octs, phits, enlts)
+                llts, nspin, iswts, octs, phits, enlts, paw_energy)
-      call elsdps_paw ( )
+     call elsdps_paw()
   endif
 
   if (file_recon.ne.' ')  call write_paw_recon ( )

atomic/set_rho_core.f90

@@ -140,7 +140,7 @@ subroutine set_rho_core
            psccharge(1:grid%mesh) = 0._dp
         end if
      else
-        aeccharge(1:grid%mesh) = rhoc(1:grid%mesh)
+!        aeccharge(1:grid%mesh) = rhoc(1:grid%mesh)
         if (nlcc) then
            psccharge(1:grid%mesh) = rhoc(1:grid%mesh)
         else