diff --git a/PW/setup.f90 b/PW/setup.f90
index b5a065942..5d9ec0583 100644
--- a/PW/setup.f90
+++ b/PW/setup.f90
@@ -97,7 +97,8 @@ SUBROUTINE setup()
       tipo,           &!
       is,             &!
       nb,             &!
-      nbe,            &
+      nbe,            &!
+      ind, ind1,      &!
       l,              &!
       ibnd             !
   LOGICAL :: &
@@ -431,15 +432,26 @@ SUBROUTINE setup()
               !
               IF ( l /= 0 ) THEN
                  !
-                 vionl = ( ( l + 1.D0 ) * dion(nbe+1,nbe+1,nt) + &
-                      l * dion(nbe,nbe,nt) ) / ( 2.D0 * l + 1.D0 )
+                 IF (ABS(jjj(nbe,nt)-lll(nbe,nt)+0.5d0).LT.1.d-7) THEN
+                    IF (ABS(jjj(nbe+1,nt)-lll(nbe+1,nt)-0.5d0).GT.1.d-7) &
+                       call errore('setup','wrong beta functions',1)
+                    ind=nbe+1
+                    ind1=nbe
+                 ELSE
+                    IF (ABS(jjj(nbe+1,nt)-lll(nbe+1,nt)+0.5d0).GT.1.d-7) &
+                       call errore('setup','wrong beta functions',1)
+                    ind=nbe
+                    ind1=nbe+1
+                 ENDIF
+                 !
+                 vionl = ( ( l + 1.D0 ) * dion(ind,ind,nt) + &
+                      l * dion(ind1,ind1,nt) ) / ( 2.D0 * l + 1.D0 )
                  !
                  betar(1:mesh(nt),nb,nt) = 1.D0 / ( 2.D0 * l + 1.D0 ) * &
-                      ( ( l + 1.D0 ) * SQRT( dion(nbe+1,nbe+1,nt) / vionl ) * &
-                      betar(1:mesh(nt),nbe+1,nt) + &
-                      l * SQRT( dion(nbe,nbe,nt) / vionl ) * &
-                      betar(1:mesh(nt),nbe,nt) )
-
+                      ( ( l + 1.D0 ) * SQRT( dion(ind,ind,nt) / vionl ) * &
+                      betar(1:mesh(nt),ind,nt) + &
+                      l * SQRT( dion(ind1,ind1,nt) / vionl ) * &
+                      betar(1:mesh(nt),ind1,nt) )
                  !
                  dion(nb,nb,nt) = vionl
                  !
@@ -453,6 +465,8 @@ SUBROUTINE setup()
                  !
               END IF
               !
+              lll(nb,nt)=lll(nbe,nt)
+              !
            END DO
            !
            nbe = 0
@@ -462,7 +476,7 @@ SUBROUTINE setup()
               nbe = nbe + 1
               !
               IF ( lchi(nb,nt) /= 0 .AND. &
-                   ABS( jchi(nb,nt) - lchi(nb,nt) - 0.5D0 ) < 1.D-7 ) nbe = nbe - 1
+                   ABS(jchi(nb,nt)-lchi(nb,nt)-0.5D0 ) < 1.D-7 ) nbe = nbe - 1
               !
            END DO
            !
@@ -478,8 +492,20 @@ SUBROUTINE setup()
               !
               IF ( l /= 0 ) THEN
                  !
-                 chi(1:mesh(nt),nb,nt)=( ( l + 1.D0 ) * chi(1:mesh(nt),nbe+1,nt)+ &
-                      l * chi(1:mesh(nt),nbe,nt)) / ( 2.D0 * l + 1.D0 )
+                 IF (ABS(jchi(nbe,nt)-lchi(nbe,nt)+0.5d0).LT.1.d-7) THEN
+                    IF (ABS(jchi(nbe+1,nt)-lchi(nbe+1,nt)-0.5d0).GT.1.d-7) &
+                       call errore('setup','wrong chi functions',1)
+                    ind=nbe+1
+                    ind1=nbe
+                 ELSE
+                    IF (ABS(jchi(nbe+1,nt)-lchi(nbe+1,nt)+0.5d0).GT.1.d-7) &
+                       call errore('setup','wrong chi functions',1)
+                    ind=nbe
+                    ind1=nbe+1
+                 END IF
+                 !
+                 chi(1:mesh(nt),nb,nt)=((l+1.D0) * chi(1:mesh(nt),ind,nt)+ &
+                      l * chi(1:mesh(nt),ind1,nt)) / ( 2.D0 * l + 1.D0 )
 
                  nbe = nbe + 1
                  !
@@ -489,6 +515,8 @@ SUBROUTINE setup()
                  !
               END IF
               !
+              lchi(nb,nt)= lchi(nbe,nt)
+              !
            END DO
            !
         END IF