diff --git a/CPV/src/Makefile b/CPV/src/Makefile
index cab341c75042797b27ac3ef1052e07679b458fb7..18e66c8dfeebfafad2de4817aa5a4cd400230d2d 100644
--- a/CPV/src/Makefile
+++ b/CPV/src/Makefile
@@ -5,6 +5,7 @@ include	../../make.inc
 # location of needed modules and included files (if any)
 MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG)../../KS_Solvers/Davidson
 
 FOBJS = \
@@ -119,7 +120,7 @@ entropy.o
 
 QEMODS=../../Modules/libqemod.a
 QEFFT=../../FFTXlib/libqefft.a
-QELA=../../KS_Solvers/CG/libcg.a ../../KS_Solvers/Davidson/libdavid.a ../../LAXlib/libqela.a
+QELA=../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a ../../KS_Solvers/Davidson/libdavid.a ../../LAXlib/libqela.a
 UTIL=../../UtilXlib/libutil.a
 
 TLDEPS= bindir mods libs libiotk
diff --git a/Makefile b/Makefile
index af3b952e57d57590e141f4bd3fc97aa3cd744850..ccf509ad844a1cf69df3a5a86e90797d43a85baf 100644
--- a/Makefile
+++ b/Makefile
@@ -172,6 +172,9 @@ libdavid : libs libutil libla
 libcg : libs libutil libla
 	( cd KS_Solvers/CG ; $(MAKE) TLDEPS= all || exit 1 )
 
+libppcg : libs libutil libla
+	( cd KS_Solvers/PPCG ; $(MAKE) TLDEPS= all || exit 1 )
+
 libla : liblapack libutil libcuda
 	( cd LAXlib ; $(MAKE) TLDEPS= all || exit 1 )
 
@@ -292,7 +295,7 @@ clean :
 	touch make.inc 
 	for dir in \
 		CPV LAXlib FFTXlib UtilXlib Modules PP PW EPW \
-                KS_Solvers/CG KS_Solvers/Davidson KS_Solvers/Davidson_RCI \
+                KS_Solvers/PPCG KS_Solvers/CG KS_Solvers/Davidson KS_Solvers/Davidson_RCI \
 		NEB ACFDT COUPLE GWW XSpectra PWCOND dft-d3 \
 		atomic clib LR_Modules pwtools upftools \
 		dev-tools extlibs Environ TDDFPT PHonon GWW \
diff --git a/Modules/bcast_qes_module.f90 b/Modules/bcast_qes_module.f90
index adb03d34ee5220b2f868719eae27abceabc9830e..f5fbef0e9321b27e74fb4267701b6cace7d4a427 100644
--- a/Modules/bcast_qes_module.f90
+++ b/Modules/bcast_qes_module.f90
@@ -1002,8 +1002,11 @@ CALL mp_bcast(o%tbeta_smoothing, ionode_id, comm)
 CALL mp_bcast(o%diago_thr_init, ionode_id, comm)
 CALL mp_bcast(o%diago_full_acc, ionode_id, comm)
 CALL mp_bcast(o%diago_cg_maxiter_ispresent, ionode_id, comm)
+CALL mp_bcast(o%diago_ppcg_maxiter_ispresent, ionode_id, comm)
 IF (o%diago_cg_maxiter_ispresent) &
    CALL mp_bcast(o%diago_cg_maxiter, ionode_id, comm)
+IF (o%diago_ppcg_maxiter_ispresent) &
+   CALL mp_bcast(o%diago_ppcg_maxiter, ionode_id, comm)
 CALL mp_bcast(o%diago_david_ndim_ispresent, ionode_id, comm)
 IF (o%diago_david_ndim_ispresent) &
    CALL mp_bcast(o%diago_david_ndim, ionode_id, comm)
diff --git a/Modules/control_flags.f90 b/Modules/control_flags.f90
index 9cba551dddac288398cd3bbf85877b5a364950a3..46b05155339d9cc65b0daf8c79a572cdf908bc37 100644
--- a/Modules/control_flags.f90
+++ b/Modules/control_flags.f90
@@ -202,9 +202,10 @@ MODULE control_flags
   REAL(DP), PUBLIC  :: &
     ethr               ! the convergence threshold for eigenvalues
   INTEGER, PUBLIC :: &
-    isolve,           &! Davidson or CG or ParO diagonalization
+    isolve,           &! index selecting Davidson,  CG , PPCG or ParO diagonalization
     david,            &! max dimension of subspace in Davidson diagonalization
-    max_cg_iter        ! maximum number of iterations in a CG call
+    max_cg_iter,      &! maximum number of iterations in a CG call
+    max_ppcg_iter      ! maximum number of iterations in a PPCG call
   LOGICAL, PUBLIC :: &
     diago_full_acc = .FALSE. ! if true,  empty eigenvalues have the same
                              ! accuracy of the occupied ones
diff --git a/Modules/input_parameters.f90 b/Modules/input_parameters.f90
index eaf5bc9ade907f8728493acac8fef6cb2a12de3c..81a75e317b93cb76cba50d78a9447ce43da6ca02 100644
--- a/Modules/input_parameters.f90
+++ b/Modules/input_parameters.f90
@@ -846,7 +846,7 @@ MODULE input_parameters
           ! used only in PWscf
 
         CHARACTER(len=80) :: diagonalization = 'david'
-          ! diagonalization = 'david' or 'cg'
+          ! diagonalization = 'david', 'cg' or 'ppcg'
           ! algorithm used by PWscf for iterative diagonalization
 
         REAL(DP) :: diago_thr_init = 0.0_DP
@@ -857,6 +857,10 @@ MODULE input_parameters
           ! max number of iterations for the first iterative diagonalization
           ! using conjugate-gradient algorithm - used only in PWscf
 
+        INTEGER :: diago_ppcg_maxiter = 100
+          ! max number of iterations for the first iterative diagonalization
+          ! using projected preconditioned conjugate-gradient algorithm - used only in PWscf
+
         INTEGER :: diago_david_ndim = 4
           ! dimension of the subspace used in Davidson diagonalization
           ! used only in PWscf
diff --git a/Modules/mp_global.f90 b/Modules/mp_global.f90
index 0a60124a9e36fdbb573a4885a8d3acf81cc2f03a..3faab76037fa1ae81d4ac3270677e6d40587a231 100644
--- a/Modules/mp_global.f90
+++ b/Modules/mp_global.f90
@@ -114,6 +114,7 @@ CONTAINS
     CALL mp_start_diag ( ndiag_, world_comm, my_comm, do_distr_diag_inside_bgrp )
     !
     call set_mpi_comm_4_cg( intra_pool_comm, intra_bgrp_comm, inter_bgrp_comm )
+    call set_mpi_comm_4_ppcg( intra_pool_comm, intra_bgrp_comm, inter_bgrp_comm )
     call set_mpi_comm_4_davidson( intra_pool_comm, intra_bgrp_comm, inter_bgrp_comm )
     !
     RETURN
@@ -127,6 +128,7 @@ CONTAINS
     USE mp, ONLY : mp_comm_free
     !
     CALL unset_mpi_comm_4_cg()
+    CALL unset_mpi_comm_4_ppcg()
     CALL unset_mpi_comm_4_davidson()
     CALL mp_comm_free ( intra_bgrp_comm )
     CALL mp_comm_free ( inter_bgrp_comm )
diff --git a/Modules/qes_libs.f90 b/Modules/qes_libs.f90
index 454dde8d45e8888b593ff5bcc5aebac6f34be54c..60bb39bece55aa70acbc18744b5ffed0d3d38cae 100644
--- a/Modules/qes_libs.f90
+++ b/Modules/qes_libs.f90
@@ -3137,6 +3137,9 @@ SUBROUTINE qes_write_electron_control(xp, obj)
       CALL xml_NewElement(xp, 'diago_cg_maxiter')
          CALL xml_addCharacters(xp, obj%diago_cg_maxiter)
       CALL xml_EndElement(xp, 'diago_cg_maxiter')
+      CALL xml_NewElement(xp, 'diago_ppcg_maxiter')
+         CALL xml_addCharacters(xp, obj%diago_ppcg_maxiter)
+      CALL xml_EndElement(xp, 'diago_ppcg_maxiter')
    CALL xml_EndElement(xp, TRIM(obj%tagname))
    !
 END SUBROUTINE qes_write_electron_control
@@ -3144,7 +3147,7 @@ END SUBROUTINE qes_write_electron_control
 SUBROUTINE qes_init_electron_control(obj, tagname, diagonalization, mixing_mode, &
                               mixing_beta, conv_thr, mixing_ndim, max_nstep, real_space_q, &
                               tq_smoothing, tbeta_smoothing, diago_thr_init, &
-                              diago_full_acc, diago_cg_maxiter)
+                              diago_full_acc, diago_cg_maxiter, diago_ppcg_maxiter)
    IMPLICIT NONE
 
    TYPE(electron_control_type) :: obj
@@ -3162,6 +3165,7 @@ SUBROUTINE qes_init_electron_control(obj, tagname, diagonalization, mixing_mode,
    REAL(DP) :: diago_thr_init
    LOGICAL  :: diago_full_acc
    INTEGER  :: diago_cg_maxiter
+   INTEGER  :: diago_ppcg_maxiter
 
    obj%tagname = TRIM(tagname)
    obj%lwrite   = .TRUE.
@@ -3178,6 +3182,7 @@ SUBROUTINE qes_init_electron_control(obj, tagname, diagonalization, mixing_mode,
    obj%diago_thr_init = diago_thr_init
    obj%diago_full_acc = diago_full_acc
    obj%diago_cg_maxiter = diago_cg_maxiter
+   obj%diago_ppcg_maxiter = diago_ppcg_maxiter
 
 END SUBROUTINE qes_init_electron_control
 
diff --git a/Modules/qes_read_module.f90 b/Modules/qes_read_module.f90
index 1377a3b9667748922593372c0e575aa49b52f102..fcb6b9a293de2a6b1e4772d63b2439a520f7bbfc 100644
--- a/Modules/qes_read_module.f90
+++ b/Modules/qes_read_module.f90
@@ -4722,6 +4722,34 @@ MODULE qes_read_module
        obj%diago_cg_maxiter_ispresent = .FALSE.
     END IF
     !
+    tmp_node_list => getElementsByTagname(xml_node, "diago_ppcg_maxiter")
+    tmp_node_list_size = getLength(tmp_node_list)
+    !
+    IF (tmp_node_list_size > 1) THEN
+        IF (PRESENT(ierr) ) THEN 
+           CALL infomsg("qes_read:electron_controlType","diago_ppcg_maxiter: too many occurrences")
+           ierr = ierr + 1 
+        ELSE 
+           CALL errore("qes_read:electron_controlType","diago_ppcg_maxiter: too many occurrences",10)
+        END IF
+    END IF
+    !
+    IF (tmp_node_list_size>0) THEN
+      obj%diago_ppcg_maxiter_ispresent = .TRUE.
+      tmp_node => item(tmp_node_list, 0)
+      CALL extractDataContent(tmp_node, obj%diago_ppcg_maxiter , IOSTAT = iostat_)
+      IF ( iostat_ /= 0 ) THEN
+         IF ( PRESENT (ierr ) ) THEN 
+            CALL infomsg("qes_read:electron_controlType","error reading diago_ppcg_maxiter")
+            ierr = ierr + 1
+         ELSE 
+            CALL errore ("qes_read:electron_controlType","error reading diago_ppcg_maxiter",10)
+         END IF
+      END IF
+    ELSE
+       obj%diago_ppcg_maxiter_ispresent = .FALSE.
+    END IF
+    !
     tmp_node_list => getElementsByTagname(xml_node, "diago_david_ndim")
     tmp_node_list_size = getLength(tmp_node_list)
     !
@@ -9811,4 +9839,4 @@ MODULE qes_read_module
   END SUBROUTINE qes_read_scalarQuantity
   !
   !
-END MODULE qes_read_module
\ No newline at end of file
+END MODULE qes_read_module
diff --git a/Modules/qes_types.f90 b/Modules/qes_types.f90
index 4263a457410f990d4175654bbd097218a6285675..d7315e75cce24c2bfdb2c167a76998fef9054969 100644
--- a/Modules/qes_types.f90
+++ b/Modules/qes_types.f90
@@ -628,6 +628,8 @@ MODULE qes_types_module
     LOGICAL :: diago_full_acc
     LOGICAL  :: diago_cg_maxiter_ispresent = .FALSE.
     INTEGER :: diago_cg_maxiter
+    LOGICAL  :: diago_ppcg_maxiter_ispresent = .FALSE.
+    INTEGER :: diago_ppcg_maxiter
     LOGICAL  :: diago_david_ndim_ispresent = .FALSE.
     INTEGER :: diago_david_ndim
     !
diff --git a/Modules/qexsd_input.f90 b/Modules/qexsd_input.f90
index 3ef47ec7ae951082cba0556e60e64ae1213449db..309bf0b461888a512fc9ce02edc759f9f1dc0c6d 100644
--- a/Modules/qexsd_input.f90
+++ b/Modules/qexsd_input.f90
@@ -251,8 +251,8 @@ MODULE qexsd_input
   SUBROUTINE qexsd_init_electron_control( obj,diagonalization,mixing_mode,mixing_beta,&
                                           conv_thr, mixing_ndim, max_nstep, tqr,tq_smoothing, &
                                           tbeta_smoothing, & 
-                                          diago_thr_init, diago_full_acc, diago_cg_maxiter,&
-                                          diago_david_ndim)
+                                          diago_thr_init, diago_full_acc, &
+                                          diago_cg_maxiter, diago_ppcg_maxiter, diago_david_ndim)
   !-------------------------------------------------------------------------------------------
   !
   IMPLICIT NONE
@@ -260,8 +260,8 @@ MODULE qexsd_input
   TYPE(electron_control_type)             ::  obj
   CHARACTER(LEN=*),INTENT(IN)             :: diagonalization,mixing_mode
   REAL(DP),INTENT(IN)                     :: mixing_beta, conv_thr, diago_thr_init
-  INTEGER,INTENT(IN)                      :: mixing_ndim,max_nstep,diago_cg_maxiter,&
-                                             diago_david_ndim
+  INTEGER,INTENT(IN)                      :: mixing_ndim,max_nstep, diago_cg_maxiter, &
+                                             diago_ppcg_maxiter, diago_david_ndim
   LOGICAL,INTENT(IN)                      :: diago_full_acc,tqr, tq_smoothing, tbeta_smoothing
   !
   CHARACTER(LEN=*),PARAMETER              :: TAGNAME="electron_control"
@@ -271,7 +271,8 @@ MODULE qexsd_input
                                 conv_thr=conv_thr,mixing_ndim=mixing_ndim,max_nstep=max_nstep,&
                                 tq_smoothing= tq_smoothing, tbeta_smoothing = tbeta_smoothing,& 
                                 real_space_q=tqr,diago_thr_init=diago_thr_init,& 
-                                diago_full_acc=diago_full_acc,diago_cg_maxiter=diago_cg_maxiter)
+                                diago_full_acc=diago_full_acc,diago_cg_maxiter=diago_cg_maxiter, &
+                                diago_ppcg_maxiter=diago_ppcg_maxiter)
    !
    END SUBROUTINE qexsd_init_electron_control
    !
diff --git a/Modules/read_namelists.f90 b/Modules/read_namelists.f90
index 859c0818db517aed416d49f7bd0e4d9f6ed31656..646630f2cece65731164f938cf9b6fbcaf2883fc 100644
--- a/Modules/read_namelists.f90
+++ b/Modules/read_namelists.f90
@@ -380,6 +380,7 @@ MODULE read_namelists_module
        diagonalization = 'david'
        diago_thr_init = 0.0_DP
        diago_cg_maxiter = 20
+       diago_ppcg_maxiter = 20
        diago_david_ndim = 4
        diago_full_acc = .FALSE.
        !
@@ -977,6 +978,7 @@ MODULE read_namelists_module
        CALL mp_bcast( diagonalization,      ionode_id, intra_image_comm )
        CALL mp_bcast( diago_thr_init,       ionode_id, intra_image_comm )
        CALL mp_bcast( diago_cg_maxiter,     ionode_id, intra_image_comm )
+       CALL mp_bcast( diago_ppcg_maxiter,   ionode_id, intra_image_comm )
        CALL mp_bcast( diago_david_ndim,     ionode_id, intra_image_comm )
        CALL mp_bcast( diago_full_acc,       ionode_id, intra_image_comm )
        CALL mp_bcast( sic,                  ionode_id, intra_image_comm )
diff --git a/NEB/src/Makefile b/NEB/src/Makefile
index 77fa43d039f782baecaadcb1c0adb0eaa3ef3992..2fa7d6c92789d991086acc1c4a0baa67598652bd 100644
--- a/NEB/src/Makefile
+++ b/NEB/src/Makefile
@@ -7,6 +7,7 @@ MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG)../../PW/src \
           $(MOD_FLAG)../../KS_Solvers/Davidson \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG)../../dft-d3/
 
 NEBOBJS = \
@@ -36,7 +37,8 @@ path_variables.o \
 set_defaults.o \
 stop_run_path.o 
 
-QEMODS=../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a ../../KS_Solvers/CG/libcg.a \
+QEMODS=../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
+       ../../KS_Solvers/CG/libcg.a  ../../KS_Solvers/PPCG/libppcg.a \
        ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
 PWOBJS= ../../PW/src/libpw.a
 
diff --git a/PHonon/FD/Makefile b/PHonon/FD/Makefile
index 5358e093676294f848df6ed96db35234065c0d9e..3906229255b20779be1876c4e245a39dcb00bd5b 100644
--- a/PHonon/FD/Makefile
+++ b/PHonon/FD/Makefile
@@ -7,13 +7,15 @@ MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG)../../PW/src \
           $(MOD_FLAG)../../dft-d3 \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG) ../../KS_Solvers/Davidson
 
 FDOBJS = \
         stop_pp.o 
 
 
-QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/Davidson/libdavid.a \
+QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
+         ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a ../../KS_Solvers/Davidson/libdavid.a \
          ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
 PWOBJS = ../../PW/src/libpw.a
 
diff --git a/PHonon/Gamma/Makefile b/PHonon/Gamma/Makefile
index e4183263745199767e8b5087e58a90fc30e8ce36..09f5dbe489e756583df0a904de6287aee47f8105 100644
--- a/PHonon/Gamma/Makefile
+++ b/PHonon/Gamma/Makefile
@@ -37,7 +37,8 @@ solve_ph.o \
 writedyn.o 
 
 PWOBJS = ../../PW/src/libpw.a
-QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a ../../KS_Solvers/CG/libcg.a \
+QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
+         ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a \
          ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
 
 TLDEPS= bindir mods libs pw
diff --git a/PHonon/PH/Makefile b/PHonon/PH/Makefile
index 9f4f1bc7e5a529beecbb6f05bd6381590a21225e..b564cf9c12aa597c944c8f8e69a39d8b216b9f96 100644
--- a/PHonon/PH/Makefile
+++ b/PHonon/PH/Makefile
@@ -177,7 +177,8 @@ write_eigenvectors.o
 
 LRMODS = ../../LR_Modules/liblrmod.a
 PWOBJS = ../../PW/src/libpw.a
-QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a ../../KS_Solvers/CG/libcg.a \
+QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
+         ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a \
          ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
 
 TLDEPS= bindir libs mods libdavid libcg lrmods
diff --git a/PP/src/Makefile b/PP/src/Makefile
index 6e1cb9ab70e57451e038d6ebb81d5af7acf454a8..5bf68b57423b215c86b32879b88286219c3dd7e5 100644
--- a/PP/src/Makefile
+++ b/PP/src/Makefile
@@ -7,6 +7,7 @@ MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG)../../PW/src \
           $(MOD_FLAG)../../KS_Solvers/Davidson \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG)../../dft-d3/
 
 PPOBJS = \
@@ -50,7 +51,8 @@ write_io_header.o \
 write_hamiltonians.o 
 
 PWOBJS = ../../PW/src/libpw.a
-QEOBJS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a ../../KS_Solvers/CG/libcg.a \
+QEOBJS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
+         ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a \
          ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a
 
 MODULES = $(PWOBJS) $(QEOBJS)
diff --git a/PP/src/make.depend b/PP/src/make.depend
index 217cc093445f263c383aee5fa1abafb40c3db05c..0d39ec0ae774ab2a5d486c986dbd9196b998e0a9 100644
--- a/PP/src/make.depend
+++ b/PP/src/make.depend
@@ -114,7 +114,6 @@ chdens_bspline.o : ../../Modules/cell_base.o
 chdens_bspline.o : ../../Modules/fft_base.o
 chdens_bspline.o : ../../Modules/io_global.o
 chdens_bspline.o : ../../Modules/kind.o
-chdens_bspline.o : chdens_module.o
 chdens_module.o : ../../FFTXlib/fft_interfaces.o
 chdens_module.o : ../../FFTXlib/fft_types.o
 chdens_module.o : ../../FFTXlib/scatter_mod.o
diff --git a/PW/src/Makefile b/PW/src/Makefile
index 39c56e4487b0b56cb13d8464e88054c1f0456bb2..69d5a80c4fec19149f9eddcac2fa8e077d448be2 100644
--- a/PW/src/Makefile
+++ b/PW/src/Makefile
@@ -6,6 +6,7 @@ include ../../make.inc
 MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG)../../KS_Solvers/Davidson \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG)../../dft-d3/
 PWOBJS = \
 pwscf.o 
@@ -246,7 +247,10 @@ wannier_clean.o \
 wannier_occ.o \
 wannier_enrg.o
 
-QEMODS=../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a ../../KS_Solvers/CG/libcg.a \
+QEMODS=../../Modules/libqemod.a \
+       ../../KS_Solvers/Davidson/libdavid.a \
+       ../../KS_Solvers/CG/libcg.a \
+       ../../KS_Solvers/PPCG/libppcg.a \
        ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a ../../dft-d3/libdftd3qe.a
 
 TLDEPS=bindir mods libs liblapack libblas dftd3
diff --git a/PW/src/c_bands.f90 b/PW/src/c_bands.f90
index 26e26a4511d6f09fc82b93df9281f6fba55815b7..6e79c344fa8d744b0349598f439cc4424a873541 100644
--- a/PW/src/c_bands.f90
+++ b/PW/src/c_bands.f90
@@ -66,6 +66,8 @@ SUBROUTINE c_bands( iter )
      WRITE( stdout, '(5X,"Davidson diagonalization with overlap")' )
   ELSE IF ( isolve == 1 ) THEN
      WRITE( stdout, '(5X,"CG style diagonalization")')
+  ELSE IF ( isolve == 2 ) THEN
+     WRITE( stdout, '(5X,"PPCG style diagonalization")')
   ELSE
      CALL errore ( 'c_bands', 'invalid type of diagonalization', isolve)
   END IF
@@ -141,6 +143,7 @@ SUBROUTINE diag_bands( iter, ik, avg_iter )
   ! ... Two types of iterative diagonalizations are currently used:
   ! ... a) Davidson algorithm (all-band)
   ! ... b) Conjugate Gradient (band-by-band)
+  ! ... b) Projected Preconditioned Conjugate Gradient (block)
   ! ...
   ! ... internal procedures :
   !
@@ -158,7 +161,7 @@ SUBROUTINE diag_bands( iter, ik, avg_iter )
   USE uspp,                 ONLY : vkb, nkb, okvan
   USE gvect,                ONLY : gstart
   USE wvfct,                ONLY : g2kin, nbndx, et, nbnd, npwx, btype
-  USE control_flags,        ONLY : ethr, lscf, max_cg_iter, isolve, &
+  USE control_flags,        ONLY : ethr, lscf, max_cg_iter, max_ppcg_iter, isolve, &
                                    gamma_only, use_para_diag
   USE noncollin_module,     ONLY : noncolin, npol
   USE wavefunctions, ONLY : evc
@@ -179,7 +182,7 @@ SUBROUTINE diag_bands( iter, ik, avg_iter )
   !
   REAL (KIND=DP), INTENT(INOUT) :: avg_iter
   !
-  REAL (KIND=DP) :: cg_iter
+  REAL (KIND=DP) :: cg_iter, ppcg_iter
   ! (weighted) number of iterations in Conjugate-Gradient
   INTEGER :: npw, ig, dav_iter, ntry, notconv
   ! number of iterations in Davidson
@@ -190,6 +193,8 @@ SUBROUTINE diag_bands( iter, ik, avg_iter )
   LOGICAL :: lrot
   ! .TRUE. if the wfc have already be rotated
   !
+  integer, parameter :: sbsize = 5, rrstep = 7 ! block dimensions used in PPCG 
+  !
 ! Davidson diagonalization uses these external routines on groups of nvec bands
   external h_psi, s_psi, g_psi
 ! subroutine h_psi(npwx,npw,nvec,psi,hpsi)  computes H*psi
@@ -201,6 +206,10 @@ SUBROUTINE diag_bands( iter, ik, avg_iter )
 !  subroutine h_1psi(npwx,npw,psi,hpsi,spsi)  computes H*psi and S*psi
 !  subroutine s_1psi(npwx,npw,psi,spsi)  computes S*psi (if needed)
 ! In addition to the above ithe initial wfc rotation uses h_psi, and s_psi
+!------------------------------------------------------------------------
+! PPCG diagonalization uses these external routines on groups of bands
+! subroutine h_psi(npwx,npw,nvec,psi,hpsi)  computes H*psi
+! subroutine s_psi(npwx,npw,nvec,psi,spsi)  computes S*psi (if needed)
 
   ALLOCATE( h_diag( npwx, npol ), STAT=ierr )
   IF( ierr /= 0 ) &
@@ -262,9 +271,9 @@ CONTAINS
     !
     IMPLICIT NONE
     !
-    IF ( isolve == 1 ) THEN
+    IF ( isolve == 1 .OR. isolve == 2 ) THEN
        !
-       ! ... Conjugate-Gradient diagonalization
+       ! ... (Projected Preconditioned) Conjugate-Gradient diagonalization
        !
        ! ... h_diag is the precondition matrix
        !
@@ -288,11 +297,22 @@ CONTAINS
              !
           END IF
           !
-          CALL rcgdiagg( h_1psi, s_1psi, h_diag, &
+          IF ( isolve == 1 ) THEN
+             CALL rcgdiagg( h_1psi, s_1psi, h_diag, &
                          npwx, npw, nbnd, evc, et(1,ik), btype(1,ik), &
                          ethr, max_cg_iter, .NOT. lscf, notconv, cg_iter )
+             !
+             avg_iter = avg_iter + cg_iter
+             !
+          ELSE
+             CALL ppcg_gamma( h_psi, s_psi, okvan, h_diag, &
+                         npwx, npw, nbnd, evc, et(1,ik), btype(1,ik), &
+                         ethr, max_ppcg_iter, notconv, ppcg_iter, sbsize , rrstep, iter  )
+             !
+             avg_iter = avg_iter + cg_iter
+             !
+          END IF
           !
-          avg_iter = avg_iter + cg_iter
           !
           ntry = ntry + 1
           !
@@ -332,7 +352,7 @@ CONTAINS
              CALL regterg (  h_psi, s_psi, okvan, g_psi, &
                          npw, npwx, nbnd, nbndx, evc, ethr, &
                          et(1,ik), btype(1,ik), notconv, lrot, dav_iter ) !    BEWARE gstart has been removed from call
-          END IF
+          ENDIF
           !
           avg_iter = avg_iter + dav_iter
           !
@@ -396,10 +416,10 @@ CONTAINS
        !
     END IF
     !
-    !write (*,*) ' current isolve value ( 1 CG, 2 Davidson)', isolve; FLUSH(6)
-    IF ( isolve == 1 ) THEN
+    !write (*,*) ' current isolve value ( 0 Davidson, 1 CG, 2 PPCG)', isolve; FLUSH(6)
+    IF ( isolve == 1 .OR. isolve == 2) THEN
        !
-       ! ... Conjugate-Gradient diagonalization
+       ! ... (Projected Preconditioned) Conjugate-Gradient diagonalization
        !
        ! ... h_diag is the precondition matrix
        !
@@ -426,12 +446,22 @@ CONTAINS
              !
           END IF
           !
-          CALL ccgdiagg( h_1psi, s_1psi, h_diag, &
+          IF ( isolve == 1) then
+             CALL ccgdiagg( h_1psi, s_1psi, h_diag, &
                          npwx, npw, nbnd, npol, evc, et(1,ik), btype(1,ik), &
                          ethr, max_cg_iter, .NOT. lscf, notconv, cg_iter )
-          !
-          avg_iter = avg_iter + cg_iter
-          !
+             !
+             avg_iter = avg_iter + cg_iter
+             !
+          ELSE
+! BEWARE npol should be added to the arguments
+             CALL ppcg_k( h_psi, s_psi, okvan, h_diag, &
+                         npwx, npw, nbnd, evc, et(1,ik), btype(1,ik), &
+                         ethr, max_ppcg_iter, notconv, ppcg_iter, sbsize , rrstep, iter)
+             !
+             avg_iter = avg_iter + ppcg_iter
+             !
+          END IF
           ntry = ntry + 1
           !
           ! ... exit condition
diff --git a/PW/src/init_run.f90 b/PW/src/init_run.f90
index 651f77fdb65b2e189c4095a96afb8dbe7fec1ac3..e4ff9001ed76353ff0cfa14238d1b4886bfd0215 100644
--- a/PW/src/init_run.f90
+++ b/PW/src/init_run.f90
@@ -67,7 +67,7 @@ SUBROUTINE init_run()
   CALL ggens( dffts, gamma_only, at, g, gg, mill, gcutms, ngms )
   if (gamma_only) THEN
      ! ... Solvers need to know gstart
-     call export_gstart_2_cg(gstart); call export_gstart_2_davidson(gstart)
+     call export_gstart_2_cg(gstart); call export_gstart_2_ppcg(gstart); call export_gstart_2_davidson(gstart)
   END IF
   !
   IF (do_comp_esm) CALL esm_init()
diff --git a/PW/src/input.f90 b/PW/src/input.f90
index 335206418c2927760cfa2b6ff1cf5145a6a5ed86..f8bad00089b01e53409f8098b33e6c8de49080c8 100644
--- a/PW/src/input.f90
+++ b/PW/src/input.f90
@@ -147,7 +147,7 @@ SUBROUTINE iosys()
   USE relax,         ONLY : epse, epsf, epsp, starting_scf_threshold
   !
   USE extrapolation, ONLY : pot_order, wfc_order
-  USE control_flags, ONLY : isolve, max_cg_iter, david, tr2, imix, gamma_only,&
+  USE control_flags, ONLY : isolve, max_cg_iter, max_ppcg_iter, david, tr2, imix, gamma_only,&
                             nmix, iverbosity, smallmem, niter, &
                             io_level, ethr, lscf, lbfgs, lmd, &
                             lbands, lconstrain, restart, twfcollect, &
@@ -264,7 +264,8 @@ SUBROUTINE iosys()
   USE input_parameters, ONLY : electron_maxstep, mixing_mode, mixing_beta, &
                                mixing_ndim, mixing_fixed_ns, conv_thr,     &
                                tqr, tq_smoothing, tbeta_smoothing,         &
-                               diago_thr_init, diago_cg_maxiter,           &
+                               diago_thr_init,                             &
+                               diago_cg_maxiter, diago_ppcg_maxiter,       &
                                diago_david_ndim, diagonalization,          &
                                diago_full_acc, startingwfc, startingpot,   &
                                real_space, scf_must_converge
@@ -932,15 +933,20 @@ SUBROUTINE iosys()
   ENDIF
   !
   SELECT CASE( trim( diagonalization ) )
+  CASE ( 'david', 'davidson' )
+     !
+     isolve = 0
+     david = diago_david_ndim
+     !
   CASE ( 'cg' )
      !
      isolve = 1
      max_cg_iter = diago_cg_maxiter
      !
-  CASE ( 'david', 'davidson' )
+  CASE ( 'ppcg' )
      !
-     isolve = 0
-     david = diago_david_ndim
+     isolve = 2
+     max_ppcg_iter = diago_ppcg_maxiter
      !
   CASE DEFAULT
      !
diff --git a/PW/src/pw_init_qexsd_input.f90 b/PW/src/pw_init_qexsd_input.f90
index eac699687ea9b85b1b2373b5022e17328fccbe5b..161dc955715faf8c25314b3f221db2a928022307 100644
--- a/PW/src/pw_init_qexsd_input.f90
+++ b/PW/src/pw_init_qexsd_input.f90
@@ -37,7 +37,8 @@
                                 ip_nr2b=>nr2b, ip_nr3b => nr3b,                                                       &
                                 ip_diagonalization=>diagonalization, mixing_mode, mixing_beta,                        &
                                 mixing_ndim, tqr, tq_smoothing, tbeta_smoothing, electron_maxstep,                    &
-                                diago_thr_init, diago_full_acc, diago_cg_maxiter, diago_david_ndim,                   &
+                                diago_thr_init, diago_full_acc,                                                       & 
+                                diago_cg_maxiter, diago_ppcg_maxiter, diago_david_ndim,                               &
                                 nk1, nk2, nk3, k1, k2, k3, nkstot, ip_xk => xk, ip_wk => wk,                          &
                                 ion_dynamics, upscale, remove_rigid_rot, refold_pos, pot_extrapolation,               &
                                 wfc_extrapolation, ion_temperature, tempw, tolp, delta_t, nraise, ip_dt => dt,        &
@@ -247,7 +248,7 @@
   END IF
   CALL qexsd_init_electron_control(obj%electron_control, diagonalization, mixing_mode, mixing_beta, conv_thr,         &
                                    mixing_ndim, electron_maxstep, tqr, tq_smoothing, tbeta_smoothing, diago_thr_init, & 
-                                   diago_full_acc, diago_cg_maxiter,  diago_david_ndim )
+                                   diago_full_acc, diago_cg_maxiter,  diago_ppcg_maxiter, diago_david_ndim )
   !--------------------------------------------------------------------------------------------------------------------------------
   !                                                   K POINTS IBZ ELEMENT
   !------------------------------------------------------------------------------------------------------------------------------ 
diff --git a/PW/tools/Makefile b/PW/tools/Makefile
index 1196671f44090ddcfe91e5fe982506b8b21dd809..9b75fb25b4d0d3cb2130bd0f51cbae4dc85c5182 100644
--- a/PW/tools/Makefile
+++ b/PW/tools/Makefile
@@ -6,8 +6,10 @@ include ../../make.inc
 MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG)../../KS_Solvers/Davidson \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG)../src
-QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/Davidson/libdavid.a \
+QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a \
+         ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a ../../KS_Solvers/Davidson/libdavid.a \
          ../../LAXlib/libqela.a ../../UtilXlib/libutil.a
 PWOBJS = ../src/libpw.a
 
diff --git a/PWCOND/src/Makefile b/PWCOND/src/Makefile
index 95cb0f1ebf43887f18d35f99db25ccbba82fb345..92c4f4c7f2911b3f1816a800d3d61d360ca5c181 100644
--- a/PWCOND/src/Makefile
+++ b/PWCOND/src/Makefile
@@ -6,6 +6,7 @@ include ../../make.inc
 MODFLAGS= $(BASEMOD_FLAGS) \
           $(MOD_FLAG) ../../KS_Solvers/Davidson \
           $(MOD_FLAG)../../KS_Solvers/CG \
+          $(MOD_FLAG)../../KS_Solvers/PPCG \
           $(MOD_FLAG)../../PW/src \
           $(MOD_FLAG)../../dft-d3
 
@@ -49,7 +50,8 @@ sunitary.o \
 transmit.o 
 
 PWOBJS = ../../PW/src/libpw.a
-QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/Davidson/libdavid.a \
+QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a \
+         ../../KS_Solvers/CG/libcg.a ../../KS_Solvers/PPCG/libppcg.a \
          ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a ../../UtilXlib/libutil.a \
 	 ../../dft-d3/libdftd3qe.a
 
diff --git a/README.md b/README.md
index fd870aa82e8f62d18b07240800c10c395923b68b..1d8760059aabb19018cb8a887545f2dd9a2c3eea 100644
--- a/README.md
+++ b/README.md
@@ -36,7 +36,7 @@ The following libraries have been isolated and partially encapsulated in view of
 - UtilXlib: performing basic MPI handling, error handling, timing handling.
 - FFTXlib: parallel (MPI and OpenMP) distributed three-dimensional FFTs, performing also load-balanced distribution of data (plane waves, G-vectors and real-space grids) across processors.
 - LAXlib: parallel distributed dense-matrix diagonalization, using ELPA, SCALapack, or a custom algorithm.
-- KS Solver: parallel iterative diagonalization for the Kohn-Sham Hamiltonian (represented as an operator),using block Davidson and band-by-band Conjugate-Gradient algorithms.
+- KS Solver: parallel iterative diagonalization for the Kohn-Sham Hamiltonian (represented as an operator),using block Davidson and band-by-band or block Conjugate-Gradient algorithms.
 - LRlib: performs a variety of tasks connected with (time-dependent) DFPT, to be used also in connection with Many-Body Perturbation Theory.
 
 ## Contributing
diff --git a/TDDFPT/src/make.depend b/TDDFPT/src/make.depend
index 22d103347cc0599bba716bcb72ea026929311cd6..6bb61ac66352df99dfe68ee2a2d3a8763f8109ed 100644
--- a/TDDFPT/src/make.depend
+++ b/TDDFPT/src/make.depend
@@ -453,6 +453,7 @@ lr_run_nscf.o : ../../Modules/cell_base.o
 lr_run_nscf.o : ../../Modules/control_flags.o
 lr_run_nscf.o : ../../Modules/fft_base.o
 lr_run_nscf.o : ../../Modules/io_files.o
+lr_run_nscf.o : ../../Modules/io_global.o
 lr_run_nscf.o : ../../Modules/mp_bands.o
 lr_run_nscf.o : ../../Modules/recvec.o
 lr_run_nscf.o : ../../PW/src/atomic_wfc_mod.o
diff --git a/atomic/src/make.depend b/atomic/src/make.depend
index 91d5af4a4bdf9736a7d50abaa7910dc70960b2c6..768d873e0e24d3a258abf8081ad79de1912358e1 100644
--- a/atomic/src/make.depend
+++ b/atomic/src/make.depend
@@ -185,6 +185,10 @@ intref.o : ../../Modules/kind.o
 intref.o : ../../Modules/radial_grids.o
 inward.o : ../../Modules/kind.o
 inward.o : ../../Modules/radial_grids.o
+kli.o : ../../Modules/kind.o
+kli.o : ../../Modules/radial_grids.o
+kli.o : ld1inc.o
+kli.o : parameters.o
 ld1.o : ../../Modules/command_line_options.o
 ld1.o : ../../Modules/environment.o
 ld1.o : ../../Modules/mp_global.o
@@ -236,6 +240,7 @@ new_potential.o : ../../Modules/constants.o
 new_potential.o : ../../Modules/funct.o
 new_potential.o : ../../Modules/kind.o
 new_potential.o : ../../Modules/radial_grids.o
+new_potential.o : kli.o
 new_potential.o : ld1inc.o
 newd_at.o : ../../Modules/kind.o
 newd_at.o : ../../Modules/radial_grids.o