diff --git a/PW/paw_init.f90 b/PW/paw_init.f90
index 77f18ce7d..b41b76f89 100644
--- a/PW/paw_init.f90
+++ b/PW/paw_init.f90
@@ -42,9 +42,9 @@ MODULE paw_init
 
   ! Called from clean_pw
   SUBROUTINE deallocate_paw_internals
-    USE paw_variables
     USE uspp_param, ONLY : upf
     USE ions_base,  ONLY : nat, ntyp => nsp
+    USE paw_variables
     !
     IMPLICIT NONE
     INTEGER :: nt, na
@@ -62,6 +62,9 @@ MODULE paw_init
         ENDDO
         DEALLOCATE(rad)
     ENDIF
+    paw_is_init = .false.
+
+   RETURN
 
   END SUBROUTINE deallocate_paw_internals
 
@@ -215,7 +218,11 @@ SUBROUTINE PAW_init_onecenter()
     INTEGER, EXTERNAL :: ldim_block, gind_block
 
     IF( paw_is_init ) THEN
+<<<<<<< paw_init.f90
+        CALL errore('PAW_init_onecenter', 'Already initialized!', 1)
+=======
         CALL infomsg('PAW_init_onecenter', 'Already initialized!')
+>>>>>>> 1.22
         RETURN
     ENDIF
     !
@@ -465,12 +472,11 @@ SUBROUTINE PAW_rad_init(l, rad)
     ! Computes weights for gaussian integrals,
     ! from numerical recipes
     SUBROUTINE weights(x,w,n)
+    USE constants, ONLY : pi, eps => eps12
     implicit none
     integer :: n, i,j,m
-    real(8), parameter :: eps=1.d-14
-    real(8) :: x(n),w(n), z,z1, p1,p2,p3,pp,pi
+    real(8) :: x(n),w(n), z,z1, p1,p2,p3,pp
     
-    pi = 4._dp*atan(1._dp)
     m=(n+1)/2
     do i=1,m
         z1 = 2._dp
diff --git a/PW/paw_onecenter.f90 b/PW/paw_onecenter.f90
index 74ed32a38..d94737ab2 100644
--- a/PW/paw_onecenter.f90
+++ b/PW/paw_onecenter.f90
@@ -111,8 +111,7 @@ SUBROUTINE PAW_potential(becsum, d, energy, e_cmp)
             ! STEP: 1 [ build rho_lm (PAW_rho_lm) ]
             NULLIFY(rho_core)
             IF (i_what == AE) THEN
-               ! to pass the atom indes is dirtyer but faster and
-               ! uses less memory than to pass only a hyperslice of the array
+               ! Compute rho spherical harmonics expansion from becsum and pfunc
                CALL PAW_rho_lm(i, becsum, upf(i%t)%paw%pfunc, rho_lm)
                ! used later for xc potential:
                rho_core => upf(i%t)%paw%ae_rho_atc
@@ -120,7 +119,7 @@ SUBROUTINE PAW_potential(becsum, d, energy, e_cmp)
                sgn = +1._dp
             ELSE
                CALL PAW_rho_lm(i, becsum, upf(i%t)%paw%ptfunc, rho_lm, upf(i%t)%qfuncl)
-               !                 optional argument for pseudo part --> ^^^
+               !          optional argument for pseudo part (aug. charge) --> ^^^
                rho_core => upf(i%t)%rho_atc ! as before
                sgn = -1._dp                 ! as before
             ENDIF
@@ -521,7 +520,7 @@ END SUBROUTINE PAW_xc_potential
 !!! in order to support non-spherical charges (as Y_lm expansion)
 !!! Note that the first derivative in vxcgc becames a gradient, while the second is a divergence.
 !!! We also have to temporary store some additional Y_lm components in order not to loose
-!!! precision during teh calculation, even if only the ones up to lmax_rho (the maximum in the
+!!! precision during the calculation, even if only the ones up to lmax_rho (the maximum in the
 !!! density of charge) matter when computing \int v * rho 
 SUBROUTINE PAW_gcxc_potential(i, rho_lm,rho_core, v_lm, energy)
     USE lsda_mod,               ONLY : nspin
@@ -869,7 +868,7 @@ SUBROUTINE PAW_h_potential(i, rho_lm, v_lm, energy)
     ! this loop computes the hartree potential using the following formula:
     !               l is the first argument in hartree subroutine
     !               r1 = min(r,r'); r2 = MAX(r,r')
-    ! V_h(r) = \sum{lm} Y_{lm}(\hat{r})/(2L+1) \int dr' 4\pi r'^2 \rho^{lm}(r') (r1^l/r2^{l+1})
+    ! V_h(r) = \sum{lm} Y_{lm}(\hat{r})/(2l+1) \int dr' 4\pi r'^2 \rho^{lm}(r') (r1^l/r2^{l+1})
     !     done here --> ^^^^^^^^^^^^^^^^^^^^^           ^^^^^^^^^^^^^^^^^^^^^^ <-- input to the hartree subroutine
     !                 output from the h.s. --> ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     DO lm = 1, i%l**2
@@ -878,7 +877,7 @@ SUBROUTINE PAW_h_potential(i, rho_lm, v_lm, energy)
             DO k = 1, i%m
                 aux(k) = pref * SUM(rho_lm(k,lm,1:nspin))
             ENDDO
-            !write(*,*) "yadda--_", l, 2*l+2, i%m
+            !
             CALL hartree(l, 2*l+2, i%m, g(i%t), aux(:), v_lm(:,lm))
     ENDDO
 
diff --git a/PW/summary.f90 b/PW/summary.f90
index 720ac8cd9..e66e8a91d 100644
--- a/PW/summary.f90
+++ b/PW/summary.f90
@@ -421,7 +421,6 @@ SUBROUTINE print_ps_info
   USE ions_base,       ONLY : ntyp => nsp
   USE atom,            ONLY : rgrid
   USE uspp_param,      ONLY : upf
-  USE paw_variables,   ONLY : rad, xlm
   USE funct,           ONLY : dft_is_gradient
 
   !
@@ -449,15 +448,9 @@ SUBROUTINE print_ps_info
      !
      WRITE( stdout, '(5x,A)') TRIM(upf(nt)%generated)
      !
-     IF(upf(nt)%tpawp) THEN
-        lmax = 2*upf(nt)%lmax_rho
-        IF ( dft_is_gradient() ) lmax = lmax + xlm
-        WRITE( stdout, '(5x, a, i4, a, i3)') &
-               "Setup to integrate on", ((lmax+1)*(lmax+2))/2, &
-               " directions: integral exact up to l =", lmax
+     IF(upf(nt)%tpawp) &
         WRITE( stdout, '(5x,a,a)') &
                "Shape of augmentation charge: ", TRIM(upf(nt)%paw%augshape)
-     ENDIF
      WRITE( stdout, '(5x,"Using radial grid of ", i4, " points, ", &
          &i2," beta functions with: ")') rgrid(nt)%mesh, upf(nt)%nbeta
      DO ib = 1, upf(nt)%nbeta