!{\src2tex{textfont=tt}}
!!****f* ABINIT/forstrnps
!! NAME
!! forstrnps
!!
!! FUNCTION
!! Compute nonlocal pseudopotential energy contribution to forces and/or stress tensor
!! as well as kinetic energy contribution to stress tensor.
!!
!! COPYRIGHT
!! Copyright (C) 1998-2014 ABINIT group (DCA, XG, GMR, AF, AR, MB, MT)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! INPUTS
!!  atindx1(natom)=index table for atoms, inverse of atindx
!!  cg(2,mcg)=wavefunctions (may be read from disk file)
!!  ecut=cut-off energy for plane wave basis sphere (Ha)
!!  ecutsm=smearing energy for plane wave kinetic energy (Ha)
!!  effmass=effective mass for electrons (1. in common case)
!!  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
!!  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
!!  istwfk(nkpt)=input option parameter that describes the storage of wfs
!!  kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis
!!  kpt(3,nkpt)=k points in reduced coordinates
!!  mband=maximum number of bands
!!  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
!!  mgfft=maximum size of 1D FFTs
!!  mkmem=number of k points treated by this node.
!!  mpi_enreg=informations about MPI parallelization
!!  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
!!  mpw= maximum number of plane waves
!!  my_natom=number of atoms treated by current processor
!!  natom=number of atoms in cell.
!!  nattyp(ntypat)= # atoms of each type.
!!  nband(nkpt)=number of bands at each k point
!!  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft
!!  nkpt=number of k points in Brillouin zone
!!  nloalg(5)=governs the choice of the algorithm for non-local operator.
!!  npwarr(nkpt)=number of planewaves in basis and boundary at each k
!!  nspinor=number of spinorial components of the wavefunctions
!!  nsppol=1 for unpolarized, 2 for spin-polarized
!!  nsym=number of elements in symmetry group
!!  ntypat=number of types of atoms
!!  occ(mband*nkpt*nsppol)=occupation numbers for each band over all k points
!!  optfor=1 if computation of forces is required
!!  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
!!  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
!!  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
!!  psps <type(pseudopotential_type)>=variables related to pseudopotentials
!!  rprimd(3,3)=dimensional primitive translations in real space (bohr)
!!  stress_needed=1 if computation of stress tensor is required
!!  symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless)
!!  use_gpu_cuda= 0 or 1 to know if we use cuda for nonlop call
!!  wtk(nkpt)=weight associated with each k point
!!  xred(3,natom)=reduced dimensionless atomic coordinates
!!  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
!!  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
!!
!! OUTPUT
!!  if (optfor==1)
!!   grnl(3*natom*optfor)=stores grads of nonlocal energy wrt atomic coordinates
!!  if (stress_needed==1)
!!   kinstr(6)=kinetic energy part of stress tensor (hartree/bohr^3)
!!   Store 6 unique components of symmetric 3x3 tensor in the order
!!   11, 22, 33, 32, 31, 21
!!   npsstr(6)=nonlocal pseudopotential energy part of stress tensor
!!    (hartree/bohr^3)
!!
!! PARENTS
!!      forstr
!!
!! CHILDREN
!!      meanvalue_g,metric,mkffnl,mkkpg,nonlop,ph1d3d,prep_bandfft_tabs
!!      prep_nonlop,stresssym,timab,xmpi_sum
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

subroutine forstrnps(atindx1,cg,ecut,ecutsm,effmass,eigen,electronpositron,&
&  grnl,istwfk,kg,kinstr,npsstr,kpt,mband,mcg,mgfft,mkmem,mpi_enreg,mpsang,&
&  mpw,my_natom,natom,nattyp,nband,ngfft,nkpt,nloalg,npwarr,nspinor,nsppol,nsym,&
&  ntypat,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,&
&  stress_needed,symrec,use_gpu_cuda,wtk,xred,ylm,ylmgr)

 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_profiling
 use m_xmpi
 use m_wffile
 use m_errors

 use m_header,           only : hdr_skip
 use m_electronpositron, only : electronpositron_type,electronpositron_calctype
 use m_pawtab,           only : pawtab_type
 use m_paw_ij,           only : paw_ij_type
 use m_pawcprj,          only : pawcprj_type

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'forstrnps'
 use interfaces_18_timing
 use interfaces_32_util
 use interfaces_41_geometry
 use interfaces_53_spacepar
 use interfaces_65_nonlocal
 use interfaces_66_wfs
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: mband,mcg,mgfft,mkmem,mpsang,mpw,my_natom,natom,nkpt
 integer,intent(in) :: nsppol,nspinor,nsym,ntypat,optfor,stress_needed
 integer,intent(in) :: use_gpu_cuda
 real(dp),intent(in) :: ecut,ecutsm,effmass
 type(electronpositron_type),pointer :: electronpositron
 type(MPI_type),intent(inout) :: mpi_enreg
 type(pseudopotential_type),intent(in) :: psps
!arrays
 integer,intent(in) :: atindx1(natom),istwfk(nkpt)
 integer,intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat),nband(nkpt*nsppol)
 integer,intent(in) :: ngfft(18),nloalg(5),npwarr(nkpt)
 integer,intent(in) :: symrec(3,3,nsym)
 real(dp),intent(in) :: cg(2,mcg)
 real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kpt(3,nkpt)
 real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom)
 real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom)
 real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm)
 real(dp),intent(out) :: grnl(3*natom*optfor),kinstr(6),npsstr(6)
 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)

!Local variables-------------------------------
!scalars
 integer,parameter :: tim_rwwf=7
 integer :: bdtot_index,choice,cplex_dij,cpopt,dimdij,dimekb1,dimekb2,dimffnl
 integer :: ia,iatom,iatom_tot,iband,icg,ider,idir,ierr,ii,ikg,ikpt,ilm,ilmn
 integer :: ipositron,ipw,ishift,isp,ispden,isppol,istwf_k,itypat
 integer :: matblk,me_distrb,my_nspinor,n1,n2,n3,nband_k,nkpg
 integer :: nnlout,npw_k,paw_opt,signs,spaceComm,tim_nonlop,tim_nonlop_prep
 integer :: blocksize,iblock,iblocksize,ibs,nblockbd
 real(dp) :: ar,arg,renorm_factor,dfsm,ecutsm_inv,eig_k,fact_kin,fsm,htpisq,kgc1
 real(dp) :: kgc2,kgc3,kin,ucvol,xx
 !character(len=500) :: message
!arrays
 integer :: gbound_dum(1,1)
 integer,allocatable :: kg_k(:,:)
 real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3)
 real(dp) :: nonlop_dum(1,1),rmet(3,3),tsec(2)
 real(dp) :: kinpw_dum(1,1)
 real(dp),allocatable :: cwavef(:,:),ekb(:,:,:)
 real(dp),allocatable :: enlout(:),ffnl(:,:,:,:),kpg_k(:,:),kstr1(:),kstr2(:)
 real(dp),allocatable :: kstr3(:),kstr4(:),kstr5(:),kstr6(:)
 real(dp),allocatable :: ph3d(:,:,:),phkxred(:,:),sij(:,:),ylm_k(:,:)
 real(dp),allocatable :: ylmgr_k(:,:,:)
 real(dp),allocatable :: lambda_loc(:),occblock(:),weight(:)
 type(pawcprj_type) :: cprj_dum(1,1)

!*************************************************************************

 call timab(920,1,tsec)
 call timab(921,1,tsec)

!Init mpicomm and me
 if(mpi_enreg%paral_kgb==1)then
   spaceComm=mpi_enreg%comm_kpt
   me_distrb=mpi_enreg%me_kpt
 else
!* In case of HF calculation
   if (mpi_enreg%paral_hf==1) then
     spaceComm=mpi_enreg%comm_kpt
     me_distrb=mpi_enreg%me_kpt
   else
     spaceComm=mpi_enreg%comm_cell
     me_distrb=mpi_enreg%me_cell
   end if
 end if

!Some constants
 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
 ipositron=abs(electronpositron_calctype(electronpositron))
 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
!Smearing of plane wave kinetic energy
 ecutsm_inv=zero;if( ecutsm>1.0d-20) ecutsm_inv=1/ecutsm
!htpisq is (1/2) (2 Pi) **2:
 htpisq=0.5_dp*(two_pi)**2

!Arrays initializations
 ABI_ALLOCATE(kg_k,(3,mpw))
 ABI_ALLOCATE(phkxred,(2,natom))
 grnl(:)=zero
 !if (optfor==1) grnl(:)=zero
 if (stress_needed==1) then
   kinstr(:)=zero;npsstr(:)=zero
 end if

!Compute different geometric tensor, as well as ucvol, from rprimd
 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)

!Common data for "nonlop" routine
 signs=1 ; idir=0  ; ishift=0 ; tim_nonlop=4 ; tim_nonlop_prep=12
 choice=2*optfor;if (stress_needed==1) choice=10*choice+3
 if (optfor==1.and.stress_needed==1)  ishift=6
 nnlout=max(1,6*stress_needed+3*natom*optfor)
 if (psps%usepaw==0) then
   paw_opt=0 ; cpopt=-1
 else
   paw_opt=2 ; cpopt=-1
 end if

!Non-local factors:
!Norm-conserving: kleimann-Bylander energies
!PAW: Dij coefficients and overlap coefficients
 if (psps%usepaw==0) then
   dimekb1=psps%dimekb;dimekb2=ntypat
   ABI_ALLOCATE(ekb,(dimekb1,dimekb2,nspinor**2))
   ekb(1:dimekb1,1:dimekb2,1)=psps%ekb(1:dimekb1,1:dimekb2)
   if (nspinor==2) then
     ekb(1:dimekb1,1:dimekb2,2)=psps%ekb(1:dimekb1,1:dimekb2)
     ekb(:,:,3:4)=zero
   end if
   if (ipositron==1) ekb(:,:,:)=-ekb(:,:,:)
   ABI_ALLOCATE(sij,(0,0))
 else
   cplex_dij=nspinor
   dimekb1=psps%dimekb*cplex_dij;dimekb2=natom
   ABI_ALLOCATE(ekb,(dimekb1,dimekb2,nspinor**2))
   ABI_ALLOCATE(sij,(dimekb1,ntypat))
   do itypat=1,ntypat
     if (cplex_dij==1) then
       sij(1:pawtab(itypat)%lmn2_size,itypat)=pawtab(itypat)%sij(:)
     else
       do ilmn=1,pawtab(itypat)%lmn2_size
         sij(2*ilmn-1,itypat)=pawtab(itypat)%sij(ilmn)
         sij(2*ilmn  ,itypat)=zero
       end do
     end if
   end do
 end if

 call timab(921,2,tsec)

!LOOP OVER SPINS
 bdtot_index=0;icg=0
 do isppol=1,nsppol

   call timab(927,1,tsec)

   
!  PAW: retrieve Dij coefficients for this spin component
   if (psps%usepaw==1) then
     ekb(:,:,:)=zero
     do ispden=1,nspinor**2
       do iatom=1,my_natom
         iatom_tot=iatom;if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%my_atmtab(iatom)
         isp=isppol;if (nspinor==2) isp=ispden
         dimdij=paw_ij(iatom)%cplex_dij*paw_ij(iatom)%lmn2_size
         if (dimdij>dimekb1) then
           MSG_BUG('size of paw_ij%dij>dimekb1!')
         end if
         do ilmn=1,dimdij
           ekb(ilmn,iatom_tot,ispden)=paw_ij(iatom)%dij(ilmn,isp)
         end do
       end do
     end do
     
!    In case of distribution over atomic sites
     if (mpi_enreg%nproc_atom>1) then
       call xmpi_sum(ekb,mpi_enreg%comm_atom,ierr)
     end if

   end if ! end PAW Dij retrieval

   call timab(927,2,tsec)

   !ekb = one

!  Loop over k points
   ikg=0
   do ikpt=1,nkpt

     nband_k=nband(ikpt+(isppol-1)*nkpt)
     istwf_k=istwfk(ikpt)
     npw_k=npwarr(ikpt)

     if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
       bdtot_index=bdtot_index+nband_k
       cycle
     end if

     call timab(922,1,tsec)

!    Define several block values for band parallelization
     if (mpi_enreg%paral_kgb/=1) then
       nblockbd=nband_k
       blocksize=1
     else
       nblockbd=nband_k/mpi_enreg%nproc_fft
       if (nband_k/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1
       nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
       blocksize=nband_k/nblockbd
     end if
     ABI_ALLOCATE(enlout,(nnlout*blocksize))
     enlout(:)=zero
     ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor*blocksize))

     ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
     if (stress_needed==1) then
       if (psps%useylm==1)  then
         ABI_ALLOCATE(ylmgr_k,(npw_k,3,mpsang*mpsang*psps%useylm))
       else
         ABI_ALLOCATE(ylmgr_k,(0,0,0))
       end if
       ABI_ALLOCATE(kstr1,(npw_k))
       ABI_ALLOCATE(kstr2,(npw_k))
       ABI_ALLOCATE(kstr3,(npw_k))
       ABI_ALLOCATE(kstr4,(npw_k))
       ABI_ALLOCATE(kstr5,(npw_k))
       ABI_ALLOCATE(kstr6,(npw_k))
     else
         ABI_ALLOCATE(ylmgr_k,(0,0,0))
     end if

     kpoint(:)=kpt(:,ikpt)

     kg_k(:,:) = 0

!$OMP PARALLEL DO
     do ipw=1,npw_k
       kg_k(1,ipw)=kg(1,ipw+ikg)
       kg_k(2,ipw)=kg(2,ipw+ikg)
       kg_k(3,ipw)=kg(3,ipw+ikg)
     end do

     if (psps%useylm==1) then
!$OMP PARALLEL DO COLLAPSE(2)
       do ilm=1,mpsang*mpsang
         do ipw=1,npw_k
           ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
         end do
       end do
       if (stress_needed==1) then
!$OMP PARALLEL DO COLLAPSE(2)
         do ilm=1,mpsang*mpsang
           do ii=1,3
             do ipw=1,npw_k
               ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm)
             end do
           end do
         end do
       end if
     end if


!    Prepare kinetic contribution to stress tensor (Warning : the symmetry
!    has not been broken, like in mkkin.f or kpg3.f . It should be, in order to be coherent).
     if (stress_needed==1) then
!$OMP PARALLEL DO PRIVATE(fact_kin,ipw,kgc1,kgc2,kgc3,kin,xx,fsm,dfsm) &
!$OMP&SHARED(ecut,ecutsm,ecutsm_inv,gprimd,htpisq,kg_k,kpoint,kstr1,kstr2,kstr3,kstr4,kstr5,kstr6,npw_k)
       do ipw=1,npw_k
!        Compute Cartesian coordinates of (k+G)
         kgc1=gprimd(1,1)*(kpoint(1)+kg_k(1,ipw))+&
&         gprimd(1,2)*(kpoint(2)+kg_k(2,ipw))+&
&         gprimd(1,3)*(kpoint(3)+kg_k(3,ipw))
         kgc2=gprimd(2,1)*(kpoint(1)+kg_k(1,ipw))+&
&         gprimd(2,2)*(kpoint(2)+kg_k(2,ipw))+&
&         gprimd(2,3)*(kpoint(3)+kg_k(3,ipw))
         kgc3=gprimd(3,1)*(kpoint(1)+kg_k(1,ipw))+&
&         gprimd(3,2)*(kpoint(2)+kg_k(2,ipw))+&
&         gprimd(3,3)*(kpoint(3)+kg_k(3,ipw))
         kin=htpisq* ( kgc1**2 + kgc2**2 + kgc3**2 )
         fact_kin=1.0_dp
         if(kin>ecut-ecutsm)then
           if(kin>ecut)then
             fact_kin=0.0_dp
           else
!            See the routine mkkin.f, for the smearing procedure
             xx=(ecut-kin)*ecutsm_inv
!            This kinetic cutoff smoothing function and its xx derivatives
!            were produced with Mathematica and the fortran code has been
!            numerically checked against Mathematica.
             fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
             dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
!            d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
!            &                         (-144+45*xx))))))*fsm**3
             fact_kin=fsm+kin*(-ecutsm_inv)*dfsm
           end if
         end if
         kstr1(ipw)=fact_kin*kgc1*kgc1
         kstr2(ipw)=fact_kin*kgc2*kgc2
         kstr3(ipw)=fact_kin*kgc3*kgc3
         kstr4(ipw)=fact_kin*kgc3*kgc2
         kstr5(ipw)=fact_kin*kgc3*kgc1
         kstr6(ipw)=fact_kin*kgc2*kgc1
       end do ! ipw
     end if

!    Compute (k+G) vectors (only if useylm=1)
     nkpg=3*nloalg(5)
     ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
     if (nkpg>0) then
       call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
     end if

!    Compute nonlocal form factors ffnl at all (k+G)
!    (ider=1 computes gradients needed for stress tensor)
     ider=0;idir=0;dimffnl=1
     if (stress_needed==1) then
       ider=1;dimffnl=2+2*psps%useylm
     end if

     ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
     call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
&     gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
&     psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
&     npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
&     psps%usepaw,psps%useylm,ylm_k,ylmgr_k)

!    Allocate the arrays phkxred and ph3d, compute phkxred and eventually ph3d
     do iatom=1,natom
       ia=atindx1(iatom)
       arg=two_pi*(kpoint(1)*xred(1,ia)+kpoint(2)*xred(2,ia)+kpoint(3)*xred(3,ia))
       phkxred(1,iatom)=cos(arg)
       phkxred(2,iatom)=sin(arg)
     end do
     if(nloalg(1)<=0)then
!      Only the allocation, not the precomputation.
       matblk=nloalg(4)
       ABI_ALLOCATE(ph3d,(2,npw_k,matblk))
     else
!      Here, allocation as well as precomputation
       matblk=natom
       ABI_ALLOCATE(ph3d,(2,npw_k,matblk))
       call ph1d3d(1,natom,kg_k,matblk,natom,npw_k,n1,n2,n3,phkxred,ph1d,ph3d)
     end if

     call timab(922,2,tsec)

!    Loop over bands; accumulate forces and/or stresses
!    Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
     ABI_ALLOCATE(occblock,(blocksize))
     ABI_ALLOCATE(weight,(blocksize))
     occblock=zero
     weight=zero
     do iblock=1,nblockbd

       iband=(iblock-1)*blocksize+1
       if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/= me_distrb) cycle

!      Select occupied bands
       occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
       if( abs(maxval(occblock))>=tol8 ) then
         call timab(923,1,tsec)
         weight(:)=wtk(ikpt)*occblock(:)

!        Load contribution from n,k
         cwavef(:,1:npw_k*my_nspinor*blocksize)=&
&         cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)

         call timab(923,2,tsec)
         call timab(926,1,tsec)

         if (mpi_enreg%paral_kgb/=1) then
           if (psps%usepaw==1) eig_k=eigen(iblock+bdtot_index)

           call nonlop(atindx1,choice,cpopt,cprj_dum,dimekb1,dimekb2,dimffnl,dimffnl,ekb,&
&           enlout,ffnl,ffnl,gmet,gprimd,idir,psps%indlmn,istwf_k,kg_k,kg_k,kpg_k,kpg_k,kpoint,&
&           kpoint,eig_k,psps%lmnmax,matblk,mgfft,mpi_enreg,mpsang,psps%mpssoang,natom,&
&           nattyp,ngfft,nkpg,nkpg,nloalg,nnlout,npw_k,npw_k,my_nspinor,nspinor,ntypat,0,paw_opt,&
&           phkxred,phkxred,ph1d,ph3d,ph3d,signs,sij,nonlop_dum,tim_nonlop,&
&           ucvol,psps%useylm,cwavef,cwavef,use_gpu_cuda=use_gpu_cuda)
         else
!          Transpose the ffnl, kinpw, kpg and ph3d arrays.
           call prep_bandfft_tabs(dimffnl,ffnl,gbound_dum,ikpt,kinpw_dum,kpoint,&
&           psps%lmnmax,matblk,mgfft,mkmem,mpi_enreg,nkpg,npw_k,ntypat,0,ph3d)
           ABI_ALLOCATE(lambda_loc,(blocksize))
           if (psps%usepaw==1)lambda_loc(1:blocksize)=&
&           eigen(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
           call prep_nonlop(atindx1,choice,cpopt,cprj_dum,dimekb1,dimekb2,dimffnl,ekb,&
&           enlout,gmet,gprimd,idir,ikpt,psps%indlmn,istwf_k,&
&           kpoint,lambda_loc,psps%lmnmax,matblk,blocksize,mgfft,mpi_enreg,mpsang,psps%mpssoang,natom,&
&           nattyp,ngfft,nkpg,nloalg,nnlout,npw_k,my_nspinor,nspinor,ntypat,paw_opt,&
&           phkxred,ph1d,signs,sij,nonlop_dum,tim_nonlop_prep,&
&           ucvol,psps%useylm,cwavef,cwavef,use_gpu_cuda=use_gpu_cuda)
           ABI_DEALLOCATE(lambda_loc)
         end if

         call timab(926,2,tsec)

!        Accumulate non-local contributions from n,k
         if (optfor==1) then
           do iblocksize=1,blocksize
             ibs=nnlout*(iblocksize-1)
             grnl(1:3*natom)=grnl(1:3*natom)+weight(iblocksize)*enlout(ibs+1+ishift:ibs+3*natom+ishift)
           end do
         end if
         if (stress_needed==1) then
           do iblocksize=1,blocksize
             ibs=nnlout*(iblocksize-1)
             npsstr(1:6)=npsstr(1:6) + weight(iblocksize)*enlout(ibs+1:ibs+6)
           end do
         end if

!        Accumulate stress tensor kinetic contributions
         if (stress_needed==1) then
           call timab(924,1,tsec)
           do iblocksize=1,blocksize
             call meanvalue_g(ar,kstr1,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
             kinstr(1)=kinstr(1)+weight(iblocksize)*ar
             call meanvalue_g(ar,kstr2,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
             kinstr(2)=kinstr(2)+weight(iblocksize)*ar
             call meanvalue_g(ar,kstr3,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
             kinstr(3)=kinstr(3)+weight(iblocksize)*ar
             call meanvalue_g(ar,kstr4,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
             kinstr(4)=kinstr(4)+weight(iblocksize)*ar
             call meanvalue_g(ar,kstr5,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
             kinstr(5)=kinstr(5)+weight(iblocksize)*ar
             call meanvalue_g(ar,kstr6,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
&             cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),0)
             kinstr(6)=kinstr(6)+weight(iblocksize)*ar
           end do
           call timab(924,2,tsec)
         end if

       end if
     end do ! End of loop on block of bands

     ABI_DEALLOCATE(occblock)
     ABI_DEALLOCATE(weight)
     ABI_DEALLOCATE(enlout)
     ABI_DEALLOCATE(cwavef)

!    Incremente indexes
     bdtot_index=bdtot_index+nband_k
     if (mkmem/=0) then
       icg=icg+npw_k*my_nspinor*nband_k
       ikg=ikg+npw_k
     end if

     ABI_DEALLOCATE(ffnl)
     ABI_DEALLOCATE(kpg_k)
     ABI_DEALLOCATE(ph3d)
     ABI_DEALLOCATE(ylm_k)
     if (stress_needed==1) then
       ABI_DEALLOCATE(kstr1)
       ABI_DEALLOCATE(kstr2)
       ABI_DEALLOCATE(kstr3)
       ABI_DEALLOCATE(kstr4)
       ABI_DEALLOCATE(kstr5)
       ABI_DEALLOCATE(kstr6)
     end if
     ABI_DEALLOCATE(ylmgr_k)

   end do ! End k point loop
 end do ! End loop over spins

!Parallel case: accumulate (n,k) contributions
 if( xmpi_paral==1) then
!  Forces
   if (optfor==1) then
     call timab(65,1,tsec)
     call xmpi_sum(grnl,spaceComm,ierr)
     call timab(65,2,tsec)
   end if
!  Stresses
   if (stress_needed==1) then
     call timab(65,1,tsec)
     call xmpi_sum(kinstr,spaceComm,ierr)
     call xmpi_sum(npsstr,spaceComm,ierr)
     call timab(65,2,tsec)
   end if
 end if

 call timab(925,1,tsec)

!Deallocate temporary space
 ABI_DEALLOCATE(kg_k)
 ABI_DEALLOCATE(phkxred)
 ABI_DEALLOCATE(ekb)

 !if (psps%usepaw==1)  then
 ABI_DEALLOCATE(sij)
 !end if

!Do final normalisation of nl contribution to forces
 if (optfor==1) grnl(:)=grnl(:)/ucvol

!Do final normalizations and symmetrizations of stress tensor contributions
 if (stress_needed==1) then
   renorm_factor=-(two_pi**2)/effmass/ucvol
   kinstr(:)=kinstr(:)*renorm_factor
   if (nsym>1) then
     call stresssym(gprimd,nsym,kinstr,symrec)
     call stresssym(gprimd,nsym,npsstr,symrec)
   end if
 end if
 

 call timab(925,2,tsec)
 call timab(920,2,tsec)

end subroutine forstrnps
!!***
