TABLE OF CONTENTS


ABINIT/m_nonlop_ylm [ Modules ]

[ Top ] [ Modules ]

NAME

  m_nonlop_ylm

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_nonlop_ylm
27 
28  use defs_basis
29  use defs_abitypes
30  use m_xmpi
31  use m_abicore
32  use m_errors
33 
34  use m_geometry,    only : strconv
35  use m_kg,          only : ph1d3d, mkkpg
36  use m_pawcprj,     only : pawcprj_type
37  use m_opernla_ylm, only : opernla_ylm
38  use m_opernlb_ylm, only : opernlb_ylm
39  use m_opernlc_ylm, only : opernlc_ylm
40  use m_opernld_ylm, only : opernld_ylm
41 
42  implicit none
43 
44  private

ABINIT/nonlop_ylm [ Functions ]

[ Top ] [ Functions ]

NAME

 nonlop_ylm

FUNCTION

 * Compute application of a nonlocal operator Vnl in order to get:
    - contracted elements (energy, forces, stresses, ...), if signs=1
    - a function in reciprocal space (|out> = Vnl|in>), if signs=2
   Operator Vnl, as the following general form:
    $Vnl=sum_{R,lmn,l''m''n''} {|P_{Rlmn}> Enl^{R}_{lmn,l''m''n''} <P_{Rl''m''n''}|}$
   Operator Vnl is -- in the typical case -- the nonlocal potential.
   - With norm-conserving pseudopots, $Enl^{R}_{lmn,l''m''n''}$ is the
     Kleinmann-Bylander energy $Ekb^{R}_{ln}$.
   - In a PAW calculation, $Enl^{R}_{lmn,l''m''n''}$ are the nonlocal
     coefficients to connect projectors $D_{ij}$.
   - The |P_{Rlmn}> are the projector functions.
 * Optionnaly, in case of PAW calculation, compute:
   - Application of the overlap matrix in reciprocal space
     (<in|S|in> or (I+S)|in>).
   - Application of (Vnl-lambda.S) in reciprocal space
     (<in|Vnl-lambda.S|in> and derivatives or (Vnl-lambda.S)|in>).
 * This routine uses spherical harmonics Ylm to express Vnl.

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  choice: chooses possible output:
    choice=0 => do nothing (only compute WF projected with NL projectors)
          =1 => non-local energy contribution
          =2 => 1st derivative(s) with respect to atomic position(s)
          =3 => 1st derivative(s) with respect to strain(s)
          =23=> 1st derivative(s) with respect to atomic pos. and
                1st derivative(s) with respect to atomic pos. and strains
          =4 => 2nd derivative(s) with respect to 2 atomic pos.
          =24=> 1st derivative(s) with respect to atm. pos. and
                2nd derivative(s) with respect to 2 atomic pos.
          =5 => 1st derivative(s) with respect to k wavevector, typically
                sum_ij [ |p_i> D_ij <dp_j/dk| + |dp_i/dk> D_ij < p_j| ]
          =6 => 2nd derivative(s) with respect to 2 strains and
                mixed 2nd derivative(s) with respect to strains & atomic pos.
          =51 =>right 1st derivative(s) with respect to k wavevector, typically
                sum_ij [ |p_i> D_ij <dp_j/dk| ]
          =52 =>left 1st derivative(s) with respect to k wavevector, typically
                sum_ij [ |dp_i/dk> D_ij < p_j| ]
          =53 =>twist 1st derivative(s) with respect to k, typically
                sum_ij [ |dp_i/dk_(idir+1)> D_ij <dp_j//dk_(idir-1)|
                        -|dp_i/dk_(idir-1)> D_ij <dp_j//dk_(idir+1)|]
          =54=> mixed 2nd derivative(s) with respect to atomic pos. and left k wavevector
          =55=> mixed 2nd derivative(s) with respect to strain and right k wavevector
          =7 => apply operator $\sum_i [ |p_i> <p_i| ],
                same as overlap operator with s_ij=identity (paw_opt==3 only)
          =8 => 2nd derivatives with respect to 2 k wavevectors
          =81=> partial 2nd derivatives with respect to 2 k wavevectors,
                full derivative with respect to k1, right derivative with respect to k2,
                (derivative with respect to k of choice 51), typically
                sum_ij [ |dp_i/dk1> D_ij <dp_j/dk2| + |p_i> D_ij < d2p_j/dk1dk2| ]
    Only choices 1,2,3,23,4,5,6 are compatible with useylm=0.
    Only choices 1,2,3,5,51,52,53,7,8,81 are compatible with signs=2
  cpopt=flag defining the status of cprjin%cp(:)=<Proj_i|Cnk> scalars (see below, side effects)
  dimenl1,dimenl2=dimensions of enl (see enl)
  dimekbq=1 if enl factors do not contain a exp(-iqR) phase, 2 is they do
  dimffnlin=second dimension of ffnlin (1+number of derivatives)
  dimffnlout=second dimension of ffnlout (1+number of derivatives)
  enl(cplex_enl*dimenl1,dimenl2,nspinortot**2,dimekbq)=
  ->Norm conserving : ==== when paw_opt=0 ====
                      (Real) Kleinman-Bylander energies (hartree)
                      dimenl1=lmnmax  -  dimenl2=ntypat
                      dimekbq is 2 if Enl contains a exp(-iqR) phase, 1 otherwise
  ->PAW :             ==== when paw_opt=1, 2 or 4 ====
                      (Real or complex, hermitian) Dij coefs to connect projectors
                      dimenl1=cplex_enl*lmnmax*(lmnmax+1)/2  -  dimenl2=natom
                      These are complex numbers if cplex_enl=2
                        enl(:,:,1) contains Dij^up-up
                        enl(:,:,2) contains Dij^dn-dn
                        enl(:,:,3) contains Dij^up-dn (only if nspinor=2)
                        enl(:,:,4) contains Dij^dn-up (only if nspinor=2)
                      dimekbq is 2 if Dij contains a exp(-iqR) phase, 1 otherwise
  ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used
          for the application of the nonlocal operator to the |in> vector
  ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used
          for the application of the nonlocal operator to the |out> vector
  -----------------------------------------------------------
  gprimd(3,3)=dimensional reciprocal space primitive translations
  idir=direction of the - atom to be moved in the case (choice=2,signs=2),
                        - k point direction in the case (choice=5,signs=2S)
                          for choice 53, twisted derivative involves idir+1 and idir-1
                        - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
  istwf_k=option parameter that describes the storage of wfs
  kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector
  kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector
  kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector
  kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector
  kptin(3)=k point in terms of recip. translations, for the |in> vector
  kptout(3)=k point in terms of recip. translations, for the |out> vector
  lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
         Typically lambda is the eigenvalue (or its guess)
  lmnmax=max. number of (l,m,n) components over all types of atoms
  matblk=dimension of the arrays ph3din and ph3dout
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nattyp(ntypat)=number of atoms of each type
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpgin,nkpgout=second sizes of arrays kpgin/kpgout
  nloalg(3)=governs the choice of the algorithm for nonlocal operator
  nnlout=dimension of enlout (when signs=1 and choice>0):
         ==== if paw_opt=0, 1 or 2 ====
         choice   nnlout     |  choice   nnlout
              1   1          |      51   6 (complex)
              2   3*natom    |      52   6 (complex)
              3   6          |      53   6
              4   6*natom    |      54   9*natom
             23   6+3*natom  |      55   36 (complex)
             24   9*natom    |       6   36+18*natom
              5   3          |       8   6
                             |      81   18 (complex)
         ==== if paw_opt=3 ====
         choice   nnlout
              1   1
              2   3*natom
              5   3
             51   3
             52   3
             54   9*natom
             55   36
              7   1
              8   6
             81   9
         ==== if paw_opt=4 ====
         not available
  npwin=number of planewaves for given k point, for the |in> vector
  npwout=number of planewaves for given k point, for the |out> vector
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  nspinortot=total number of spinorial components of the wavefunctions
  ntypat=number of types of atoms in cell
  paw_opt= define the nonlocal operator concerned with:
           paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
           paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
           paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
           paw_opt=3 : PAW overlap matrix (Sij)
           paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
  phkxredin(2,natom)=phase factors exp(2 pi kptin.xred)
  phkxredout(2,natom)=phase factors exp(2 pi kptout.xred)
  ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
  ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in)
  ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out)
  signs= if 1, get contracted elements (energy, forces, stress, ...)
         if 2, applies the non-local operator to a function in reciprocal space
  sij(dimenl1,ntypat*(paw_opt/3))=overlap matrix components (only if paw_opt=2, 3 or 4)
  ucvol=unit cell volume (bohr^3)
  vectin(2,npwin*nspinor)=input cmplx wavefunction coefficients <G|in>
  [cprjin_left(natom,nspinor)]=The projected input wave function <p_nlm|in_left>
    for the left wavefunction. Data are assumed to be in memory, they are NOT recalculated here.

OUTPUT

 ==== if (signs==1) ====
 --If (paw_opt==0, 1 or 2)
    enlout(nnlout)= contribution to the non-local part of the following properties:
      if choice=1 : enlout(1)             -> the energy
      if choice=2 : enlout(3*natom)       -> 1st deriv. of energy wrt atm. pos (forces)
      if choice=3 : enlout(6)             -> 1st deriv. of energy wrt strain (stresses)
      if choice=4 : enlout(6*natom)       -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
      if choice=23: enlout(6+3*natom)     -> 1st deriv. of energy wrt atm. pos (forces) and
                                             1st deriv. of energy wrt strain (stresses)
      if choice=24: enlout(9*natom)       -> 1st deriv. of energy wrt atm. pos (forces) and
                                             2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
      if choice=5 : enlout(3)             -> 1st deriv. of energy wrt k
      if choice=51: enlout(3)             -> 1st deriv. (right) of energy wrt k
      if choice=52: enlout(3)             -> 1st deriv. (left) of energy wrt k
      if choice=53: enlout(3)             -> 1st deriv. (twist) of energy wrt k
      if choice=54: enlout(18*natom)      -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
      if choice=55: enlout(36)            -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
      if choice=6 : enlout(36+18*natom)   -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
                                             2nd deriv. of energy wrt to atm. pos and strain (internal strain)
      if choice=8 : enlout(6)             -> 2nd deriv. of energy wrt 2 k
      if choice=81: enlout(9)             -> 2nd deriv.of E: full derivative w.r.t. k1, right derivative w.r.t k2
 --If (paw_opt==3)
      if choice=1 : enlout(1)             -> contribution to <c|S|c> (note: not including <c|c>)
      if choice=2 : enlout(3*natom)       -> contribution to <c|dS/d_atm.pos|c>
      if choice=51: enlout(3)             -> contribution to <c|d(right)S/d_k|c>
      if choice=52: enlout(3)             -> contribution to <c|d(left)S/d_k|c>
      if choice=54: enlout(18*natom)      -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
      if choice=55: enlout(36)            -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
      if choice=7 : enlout(1)             -> contribution to <c|sum_i[p_i><p_i]|c>
      if choice=8 : enlout(6)             -> contribution to <c|d2S/d_k1d_k2|c>
      if choice=81: enlout(9)             -> contribution to <c|dS/d_k1[d(right)d_k2]|c>
 --If (paw_opt==4)
      not available
 ==== if (signs==2) ====
 --if (paw_opt=0, 1 or 4)
    vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
                or one of its derivatives to the input vect.:
      if (choice=1)  <G|V_nonlocal|vect_in>
      if (choice=2)  <G|dV_nonlocal/d(atm. pos)|vect_in>
      if (choice=3)  <G|dV_nonlocal/d(strain)|vect_in>
      if (choice=5)  <G|dV_nonlocal/d(k)|vect_in>
      if (choice=51) <G|d(right)V_nonlocal/d(k)|vect_in>
      if (choice=52) <G|d(left)V_nonlocal/d(k)|vect_in>
      if (choice=53) <G|d(twist)V_nonlocal/d(k)|vect_in>
      if (choice=8)  <G|d2V_nonlocal/d(k)d(k)|vect_in>
      if (choice=81) <G|d[d(right)V_nonlocal/d(k)]/d(k)|vect_in>
 --if (paw_opt=2)
    vectout(2,npwout*my_nspinor*ndat)=final vector in reciprocal space:
      if (choice=1)  <G|V_nonlocal-lamdba.(I+S)|vect_in>
      if (choice=2)  <G|d[V_nonlocal-lamdba.(I+S)]/d(atm. pos)|vect_in>
      if (choice=3)  <G|d[V_nonlocal-lamdba.(I+S)]/d(strain)|vect_in>
      if (choice=5)  <G|d[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
      if (choice=51) <G|d(right)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
      if (choice=52) <G|d(left)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
      if (choice=53) <G|d(twist)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
      if (choice=8)  <G|d2[V_nonlocal-lamdba.(I+S)]/d(k)d(k)|vect_in>
      if (choice=81) <G|d[d(right[V_nonlocal-lamdba.(I+S)]/d(k)]/d(k)|vect_in>
 --if (paw_opt=3 or 4)
    svectout(2,npwout*my_nspinor*ndat)=result of the aplication of Sij (overlap matrix)
                  or one of its derivatives to the input vect.:
      if (choice=1)  <G|I+S|vect_in>
      if (choice=2)  <G|dS/d(atm. pos)|vect_in>
      if (choice=3)  <G|dS/d(strain)|vect_in>
      if (choice=5)  <G|dS/d(k)|vect_in>
      if (choice=51) <G|d(right)S/d(k)|vect_in>
      if (choice=52) <G|d(left)S/d(k)|vect_in>
      if (choice=53) <G|d(twist)S/d(k)|vect_in>
      if (choice=3)  <G|d[V_nonlocal-lamdba.(I+S)]/d(strain)|vect_in>
      if (choice=7)  <G|sum_i[p_i><p_i]|vect_in>
      if (choice=8)  <G|d2S/d(k)d(k)|vect_in>
      if (choice=81) <G|d[d(right)S/d(k)]/d(k)|vect_in>

SIDE EFFECTS

  cprjin(natom,nspinor) <type(pawcprj_type)>=projected input wave function |in> on non-local projectors
                                            =<p_lmn|in> and derivatives
                    Treatment depends on cpopt parameter:
                     if cpopt=-1, <p_lmn|in> (and derivatives)
                                  are computed here (and not saved)
                     if cpopt= 0, <p_lmn|in> are computed here and saved
                                  derivatives are eventually computed but not saved
                     if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
                                  other derivatives are eventually computed but not saved
                     if cpopt= 2  <p_lmn|in> are already in memory;
                                  first (and 2nd) derivatives are computed here and not saved
                     if cpopt= 3  <p_lmn|in> are already in memory;
                                  first derivatives are computed here and saved
                                  other derivatives are eventually computed but not saved
                     if cpopt= 4  <p_lmn|in> and first derivatives are already in memory;
                                  other derivatives are not computed (except when choice=8 or 81)
                                  This option is not compatible with choice=4,24 or 6
                     Warning: for cpopt= 1 or 3, derivatives wrt strains do not contain
                              the contribution due to the volume change;
                              i.e. <dp_lmn/dEps|in> are incomplete.

NOTES

 This application of the nonlocal operator is programmed using a direct
 implementation of spherical harmonics (Ylm). Abinit used historically
 Legendre polynomials for the application of nonlocal operator; but the
 implementation of PAW algorithm enforced the use of Ylm.

 In the case signs=1, the array vectout is not used, nor modified
 so that the same array as vectin can be used as a dummy argument;
 the same is true for the pairs npwin-npwout, ffnlin-ffnlout,
 kgin-kgout, ph3din-ph3dout, phkredin-phkxredout).

TODO

 * Complete implementation of spin-orbit

PARENTS

      nonlop

CHILDREN

      mkkpg,opernla_ylm,opernlb_ylm,opernlc_ylm,opernld_ylm,ph1d3d,strconv
      xmpi_sum

SOURCE

 325  subroutine nonlop_ylm(atindx1,choice,cpopt,cprjin,dimenl1,dimenl2,dimekbq,dimffnlin,dimffnlout,&
 326 &                      enl,enlout,ffnlin,ffnlout,gprimd,idir,indlmn,istwf_k,&
 327 &                      kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,lmnmax,matblk,mgfft,&
 328 &                      mpi_enreg,natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,nnlout,&
 329 &                      npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,phkxredin,phkxredout,ph1d,&
 330 &                      ph3din,ph3dout,signs,sij,svectout,ucvol,vectin,vectout,cprjin_left)
 331 
 332 
 333 !This section has been created automatically by the script Abilint (TD).
 334 !Do not modify the following lines by hand.
 335 #undef ABI_FUNC
 336 #define ABI_FUNC 'nonlop_ylm'
 337 !End of the abilint section
 338 
 339  implicit none
 340 
 341 !Arguments ------------------------------------
 342 !scalars
 343  integer,intent(in) :: choice,cpopt,dimenl1,dimenl2,dimekbq,dimffnlin,dimffnlout,idir
 344  integer,intent(in) :: istwf_k,lmnmax,matblk,mgfft,natom,nkpgin,nkpgout,nnlout
 345  integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,signs
 346  real(dp),intent(in) :: lambda,ucvol
 347  type(MPI_type),intent(in) :: mpi_enreg
 348 !arrays
 349  integer,intent(in) :: atindx1(natom),kgin(3,npwin)
 350  integer,intent(in),target :: indlmn(6,lmnmax,ntypat)
 351  integer,intent(in) :: kgout(3,npwout),nattyp(ntypat),ngfft(18),nloalg(3)
 352  real(dp),intent(in) :: enl(dimenl1,dimenl2,nspinortot**2,dimekbq)
 353  real(dp),intent(in),target :: ffnlin(npwin,dimffnlin,lmnmax,ntypat)
 354  real(dp),intent(in),target :: ffnlout(npwout,dimffnlout,lmnmax,ntypat)
 355  real(dp),intent(in) :: gprimd(3,3)
 356  real(dp),intent(in),target :: kpgin(npwin,nkpgin),kpgout(npwout,nkpgout)
 357  real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom)
 358  real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom)
 359  real(dp),intent(in) :: sij(dimenl1,ntypat*((paw_opt+1)/3))
 360  real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk)
 361  real(dp),intent(inout) :: vectin(:,:)
 362  real(dp),intent(out) :: enlout(:)
 363  real(dp),intent(out) :: svectout(:,:)
 364  real(dp),intent(inout) :: vectout (:,:)
 365  type(pawcprj_type),intent(inout) :: cprjin(:,:)
 366  type(pawcprj_type),optional,intent(in) :: cprjin_left(:,:)
 367 
 368 !Local variables-------------------------------
 369 !scalars
 370  integer :: choice_a,choice_b,cplex,cplex_enl,cplex_fac,ia,ia1,ia2,ia3,ia4,ia5
 371  integer :: iatm,ic,idir1,idir2,ii,ierr,ilmn,iln,ishift,ispinor,itypat,jc,mincat,mu,mua,mub,mu0
 372  integer :: n1,n2,n3,nd2gxdt,ndgxdt,ndgxdt_stored,nd2gxdtfac,ndgxdtfac
 373  integer :: nincat,nkpgin_,nkpgout_,nlmn,nu,nua1,nua2,nub1,nub2,optder
 374  real(dp) :: enlk
 375  logical :: check,testnl
 376  character(len=500) :: message
 377 !arrays
 378  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
 379  integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
 380  integer,allocatable :: cplex_dgxdt(:),cplex_d2gxdt(:)
 381  integer,ABI_CONTIGUOUS pointer :: indlmn_typ(:,:)
 382  real(dp),allocatable :: d2gxdt(:,:,:,:,:),d2gxdtfac(:,:,:,:,:),d2gxdtfac_sij(:,:,:,:,:)
 383  real(dp),allocatable :: dgxdt(:,:,:,:,:),dgxdtfac(:,:,:,:,:),dgxdtfac_sij(:,:,:,:,:)
 384  real(dp),allocatable :: ddkk(:),fnlk(:),gmet(:,:)
 385  real(dp),allocatable :: gx(:,:,:,:),gxfac(:,:,:,:),gxfac_sij(:,:,:,:),gx_left(:,:,:,:)
 386  real(dp),allocatable :: sij_typ(:),strnlk(:)
 387  real(dp),allocatable :: work1(:),work2(:),work3(:,:),work4(:,:),work5(:,:,:),work6(:,:,:),work7(:,:,:)
 388  real(dp),ABI_CONTIGUOUS pointer :: ffnlin_typ(:,:,:),ffnlout_typ(:,:,:),kpgin_(:,:),kpgout_(:,:)
 389 
 390 ! **********************************************************************
 391 
 392  DBG_ENTER("COLL")
 393 
 394 !Check consistency of arguments
 395 !==============================================================
 396 
 397 !signs=1, almost all choices
 398  if (signs==1) then
 399    if(paw_opt<3) then
 400      check=(choice==0 .or.choice==1 .or.choice==2 .or.choice==3 .or.choice==4 .or.&
 401 &     choice==23.or.choice==24.or.choice==5 .or.choice==51.or.choice==52.or.&
 402 &     choice==53.or.choice==54.or.choice==55.or.&
 403 &     choice==6 .or.choice==8 .or.choice==81)
 404    else if (paw_opt==3) then
 405      check=(choice== 0.or.choice== 1.or.choice== 2.or.choice==3.or.choice==5.or.&
 406 &     choice==23.or.choice==51.or.choice==52.or.choice==54.or.choice==55.or.&
 407 &     choice== 8.or.choice==81)
 408    else
 409      check = .false.
 410    end if
 411    ABI_CHECK(check,'BUG: choice not compatible (for signs=1)')
 412  end if
 413 
 414 !signs=2, less choices
 415  if (signs==2) then
 416    check=(choice==0.or.choice==1.or.choice==2.or.choice==3 .or.&
 417 &   choice==5.or.choice==51.or.choice==52.or.choice==53.or.choice==54.or.&
 418 &   choice==7.or.choice==8.or.choice==81)
 419    ABI_CHECK(check,'BUG: choice not compatible (for signs=2)')
 420  end if
 421 !1<=idir<=6 is required when choice=3 and signs=2
 422  if (choice==3.and.signs==2) then
 423    check=(idir>=1.and.idir<=6)
 424    ABI_CHECK(check,'BUG: choice=3 and signs=2 requires 1<=idir<=6')
 425 !1<=idir<=9 is required when choice==8/81 and signs=2
 426  else if ((choice==8.or.choice==81.or.choice==54).and.signs==2) then
 427    check=(idir>=1.and.idir<=9)
 428    ABI_CHECK(check,'BUG: choice=8/81 and signs=2 requires 1<=idir<=9')
 429  else
 430 !  signs=2 requires 1<=idir<=3 when choice>1
 431    check=(signs/=2.or.choice<=1.or.choice==7.or.(idir>=1.and.idir<=3))
 432    ABI_CHECK(check,'BUG: signs=2 requires 1<=idir<=3')
 433  end if
 434 !check allowed values for cpopt
 435  check=(cpopt>=-1.and.cpopt<=4)
 436  ABI_CHECK(check,'bad value for cpopt')
 437  check=(cpopt/=4.or.(choice/=4.and.choice/=24.and.choice/=6))
 438  ABI_CHECK(check,'BUG: cpopt=4 not allowed for 2nd derivatives')
 439  check=(cpopt/=2.or.(choice/=8.and.choice/=81))
 440  ABI_CHECK(check,'BUG: cpopt=2 not allowed for choice=8,81, use cpopt=4 instead')
 441 !check conditions for optional arguments
 442  check=((.not.present(cprjin_left)).or.(signs==1.and.choice==1))
 443  ABI_CHECK(check,'BUG: when cprjin_left is present, must have choice=1,signs=1')
 444 !protect special case choice==7
 445  check=(choice/=7.or.paw_opt==3)
 446  ABI_CHECK(check,'BUG: when choice=7, paw_opt must be 3')
 447 !spin-orbit not yet allowed
 448  check=(maxval(indlmn(6,:,:))<=1)
 449  ABI_CHECK(check,'BUG: spin-orbit not yet allowed')
 450 
 451 !Test: size of blocks of atoms
 452  mincat=min(NLO_MINCAT,maxval(nattyp))
 453  if (nloalg(2)<=0.and.mincat>matblk) then
 454    write(message, '(a,a,a,i4,a,i4,a)' ) &
 455 &   'With nloc_mem<=0, mincat must be less than matblk.',ch10,&
 456 &   'Their value is ',mincat,' and ',matblk,'.'
 457    MSG_BUG(message)
 458  end if
 459 
 460 !Define dimensions of projected scalars
 461 !==============================================================
 462 
 463 !Define some useful variables
 464  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
 465  choice_a=merge(choice,1,choice/=7);choice_b=choice_a
 466  if (cpopt>=2) choice_a=-choice_a
 467  cplex=2;if (istwf_k>1) cplex=1 !Take into account TR-symmetry
 468  cplex_enl=1;if (paw_opt>0) cplex_enl=2*dimenl1/(lmnmax*(lmnmax+1))
 469  cplex_fac=max(cplex,dimekbq)
 470  if ((nspinortot==2.or.cplex_enl==2).and.paw_opt>0.and.choice/=7) cplex_fac=2
 471 
 472 !Define dimensions of projected scalars
 473  ndgxdt=0;ndgxdtfac=0;nd2gxdt=0;nd2gxdtfac=0
 474  if (choice==2) then
 475    if (signs==1) ndgxdt=3
 476    if (signs==2) ndgxdt=1
 477    if (signs==2) ndgxdtfac=1
 478  end if
 479  if (choice==3) then
 480    if (signs==1) ndgxdt=6
 481    if (signs==2) ndgxdt=1
 482    if (signs==2) ndgxdtfac=1
 483  end if
 484  if (choice==23) then
 485    if (signs==1) ndgxdt=9
 486  end if
 487  if (choice==4) then
 488    if(signs==1) ndgxdt=3
 489    if(signs==1) ndgxdtfac=3
 490    if(signs==1) nd2gxdt=6
 491  end if
 492  if (choice==24) then
 493    if(signs==1) ndgxdt=3
 494    if(signs==1) ndgxdtfac=3
 495    if(signs==1) nd2gxdt=6
 496  end if
 497  if (choice==5) then
 498    if(signs==1) ndgxdt=3
 499    if(signs==2) ndgxdt=1
 500    if(signs==2) ndgxdtfac=1
 501  end if
 502  if (choice==51) then
 503    if(signs==1) ndgxdt=3
 504    if(signs==2) ndgxdt=1
 505    if(signs==2) ndgxdtfac=1
 506  end if
 507  if (choice==52) then
 508    if(signs==1) ndgxdt=3
 509    if(signs==2) ndgxdt=1
 510    if(signs==2) ndgxdtfac=1
 511  end if
 512  if (choice==53) then
 513    if(signs==1) ndgxdt=3
 514    if(signs==1) ndgxdtfac=3
 515    if(signs==2) ndgxdt=2
 516    if(signs==2) ndgxdtfac=2
 517  end if
 518  if (choice==54) then
 519    if(signs==1) ndgxdt=6
 520    if(signs==1) ndgxdtfac=6
 521    if(signs==1) nd2gxdt=9
 522    if(signs==2) ndgxdt=1
 523    if(signs==2) nd2gxdt=1
 524    if(signs==2) ndgxdtfac=1
 525    if(signs==2) nd2gxdtfac=1
 526  end if
 527  if (choice==55) then
 528    if(signs==1) ndgxdt=9
 529    if(signs==1) ndgxdtfac=9
 530    if(signs==1) nd2gxdt=18
 531  end if
 532  if (choice==6) then
 533    if(signs==1) ndgxdt=9
 534    if(signs==1) ndgxdtfac=9
 535    if(signs==1) nd2gxdt=54
 536  end if
 537  if (choice==8) then
 538    if(signs==1) ndgxdt=3
 539    if(signs==1) ndgxdtfac=3
 540    if(signs==1) nd2gxdt=6
 541    if(signs==2) ndgxdt=2
 542    if(signs==2) ndgxdtfac=2
 543    if(signs==2) nd2gxdt=1
 544    if(signs==2) nd2gxdtfac=1
 545  end if
 546  if (choice==81) then
 547    if(signs==1) ndgxdt=3
 548    if(signs==1) ndgxdtfac=3
 549    if(signs==1) nd2gxdt=6
 550    if(signs==2) ndgxdt=1
 551    if(signs==2) ndgxdtfac=1
 552    if(signs==2) nd2gxdt=1
 553    if(signs==2) nd2gxdtfac=1
 554  end if
 555  ABI_CHECK(ndgxdtfac<=ndgxdt,"BUG: ndgxdtfac>ndgxdt!")
 556  optder=0;if (ndgxdtfac>0) optder=1
 557  if (nd2gxdtfac>0) optder=2
 558 
 559 !Consistency tests
 560  if (cpopt==4) then
 561    if (ndgxdt>0.and.cprjin(1,1)%ncpgr<=0) then
 562      message='cprjin%ncpgr=0 not allowed with cpopt=4 and these (choice,signs) !'
 563      MSG_BUG(message)
 564    end if
 565  end if
 566  if (cpopt==1.or.cpopt==3) then
 567    if (cprjin(1,1)%ncpgr<ndgxdt) then
 568      message='should have cprjin%ncpgr>=ndgxdt with cpopt=1 or 3 !'
 569      MSG_BUG(message)
 570    end if
 571  end if
 572 
 573 !Additional steps before calculation
 574 !==============================================================
 575 
 576 !Initialize output arrays
 577  if (signs==1) then
 578    ABI_ALLOCATE(fnlk,(3*natom))
 579    ABI_ALLOCATE(ddkk,(6))
 580    ABI_ALLOCATE(strnlk,(6))
 581    enlk=zero;fnlk=zero;ddkk=zero;strnlk=zero
 582    enlout(:)=zero
 583  end if
 584  if (signs==2) then
 585    if (paw_opt==0.or.paw_opt==1.or.paw_opt==4) vectout(:,:)=zero
 586    if (paw_opt==2.and.choice==1) vectout(:,:)=-lambda*vectin(:,:)
 587    if (paw_opt==2.and.choice> 1) vectout(:,:)=zero
 588    if (paw_opt==3.or.paw_opt==4) then
 589      if (choice==1) svectout(:,:)=vectin(:,:)
 590      if (choice> 1) svectout(:,:)=zero
 591    end if
 592  end if
 593 
 594 !Eventually re-compute (k+G) vectors (and related data)
 595  nkpgin_=0
 596  if (choice==2.or.choice==54) nkpgin_=3
 597  if (signs==1) then
 598    if (choice==4.or.choice==24) nkpgin_=9
 599    if (choice==3.or.choice==23.or.choice==6) nkpgin_=3
 600    if (choice==55) nkpgin_=3
 601  end if
 602  if (nkpgin<nkpgin_) then
 603    ABI_ALLOCATE(kpgin_,(npwin,nkpgin_))
 604    call mkkpg(kgin,kpgin_,kptin,nkpgin_,npwin)
 605  else
 606    nkpgin_ = nkpgin
 607    kpgin_  => kpgin
 608  end if
 609  nkpgout_=0
 610  if ((choice==2.or.choice==54).and.signs==2) nkpgout_=3
 611  if (nkpgout<nkpgout_) then
 612    ABI_ALLOCATE(kpgout_,(npwout,nkpgout_))
 613    call mkkpg(kgout,kpgout_,kptout,nkpgout_,npwout)
 614  else
 615    nkpgout_ = nkpgout
 616    kpgout_ => kpgout
 617  end if
 618 
 619 !Big loop on atom types
 620 !==============================================================
 621 
 622  ia1=1;iatm=0
 623  do itypat=1,ntypat
 624 
 625 !  Get atom loop indices for different types:
 626    ia2=ia1+nattyp(itypat)-1;ia5=1
 627 
 628 !  Select quantities specific to the current type of atom
 629    nlmn=count(indlmn(3,:,itypat)>0)
 630 
 631 !  Test on local part
 632    testnl=(paw_opt/=0)
 633    if (paw_opt==0) testnl=any(enl(:,:,:,:)>tol10)
 634 
 635 !  Some non-local part is to be applied for that type of atom
 636    if (testnl) then
 637 
 638 !    Store some quantities depending only of the atom type
 639      ffnlin_typ => ffnlin(:,:,:,itypat)
 640      indlmn_typ => indlmn(:,:,itypat)
 641      if (signs==2) then
 642        ffnlout_typ => ffnlout(:,:,:,itypat)
 643      end if
 644      if (paw_opt>=2) then
 645        ABI_ALLOCATE(sij_typ,(nlmn*(nlmn+1)/2))
 646        if (cplex_enl==1) then
 647          do ilmn=1,nlmn*(nlmn+1)/2
 648            sij_typ(ilmn)=sij(ilmn,itypat)
 649          end do
 650        else
 651          do ilmn=1,nlmn*(nlmn+1)/2
 652            sij_typ(ilmn)=sij(2*ilmn-1,itypat)
 653          end do
 654        end if
 655      else
 656        ABI_ALLOCATE(sij_typ,(0))
 657      end if
 658 
 659 !    Loop over atoms of the same type
 660 !    ==============================================================
 661 
 662 !    Cut the sum on different atoms in blocks, to allow memory saving.
 663 !    Inner summations on atoms will be done from ia3 to ia4.
 664 !    Note: the maximum range from ia3 to ia4 is mincat (max. increment of atoms).
 665 
 666      do ia3=ia1,ia2,mincat
 667        ia4=min(ia2,ia3+mincat-1)
 668 !      Give the increment of number of atoms in this subset.
 669        nincat=ia4-ia3+1
 670 
 671 !      Prepare the phase factors if they were not already computed
 672        if (nloalg(2)<=0) then
 673          call ph1d3d(ia3,ia4,kgin,matblk,natom,npwin,n1,n2,n3,phkxredin,ph1d,ph3din)
 674        end if
 675 
 676 !      Allocate memory for projected scalars
 677        ABI_ALLOCATE(gx,(cplex,nlmn,nincat,nspinor))
 678        ABI_ALLOCATE(dgxdt,(cplex,ndgxdt,nlmn,nincat,nspinor))
 679        ABI_ALLOCATE(d2gxdt,(cplex,nd2gxdt,nlmn,nincat,nspinor))
 680        ABI_ALLOCATE(d2gxdtfac,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor))
 681        ABI_ALLOCATE(dgxdtfac,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor))
 682        ABI_ALLOCATE(gxfac,(cplex_fac,nlmn,nincat,nspinor))
 683        gx(:,:,:,:)=zero;gxfac(:,:,:,:)=zero
 684        if (ndgxdt>0) dgxdt(:,:,:,:,:)=zero
 685        if (ndgxdtfac>0) dgxdtfac(:,:,:,:,:)=zero
 686        if (nd2gxdt>0) d2gxdt(:,:,:,:,:)=zero
 687        if (nd2gxdtfac>0) d2gxdtfac(:,:,:,:,:)=zero
 688        if (paw_opt>=3) then
 689          ABI_ALLOCATE(gxfac_sij,(cplex,nlmn,nincat,nspinor))
 690          ABI_ALLOCATE(dgxdtfac_sij,(cplex,ndgxdtfac,nlmn,nincat,nspinor))
 691          ABI_ALLOCATE(d2gxdtfac_sij,(cplex,nd2gxdtfac,nlmn,nincat,nspinor))
 692          gxfac_sij(:,:,:,:)=zero
 693          if (ndgxdtfac>0) dgxdtfac_sij(:,:,:,:,:)=zero
 694          if (nd2gxdtfac>0) d2gxdtfac_sij(:,:,:,:,:) = zero
 695        else
 696          ABI_ALLOCATE(gxfac_sij,(0,0,0,0))
 697          ABI_ALLOCATE(dgxdtfac_sij,(0,0,0,0,0))
 698          ABI_ALLOCATE(d2gxdtfac_sij,(0,0,0,0,0))
 699        end if
 700 
 701 !      When istwf_k > 1, gx derivatives can be real or pure imaginary
 702 !      cplex_dgxdt(i)  = 1 if dgxdt(1,i,:,:)  is real, 2 if it is pure imaginary
 703 !      cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary
 704        ABI_ALLOCATE(cplex_dgxdt,(ndgxdt))
 705        ABI_ALLOCATE(cplex_d2gxdt,(nd2gxdt))
 706        cplex_dgxdt(:) = 1 ; cplex_d2gxdt(:) = 1
 707        if(ndgxdt > 0) then
 708          if (choice==5.or.choice==51.or.choice==52.or.choice==53.or. &
 709 &         choice==8.or.choice==81) cplex_dgxdt(:) = 2
 710          if (choice==54.and.signs==1) cplex_dgxdt(4:6) = 2
 711          if (choice==54.and.signs==2) cplex_dgxdt(:)   = 2
 712          if (choice==55.and.signs==1) cplex_dgxdt(7:9) = 2
 713        end if
 714        if(nd2gxdt > 0) then
 715          if (choice==54) cplex_d2gxdt(:) = 2
 716          if (choice==55.and.signs==1) cplex_d2gxdt(1:18)= 2
 717        end if
 718 
 719 !      Compute projection of current wave function |c> on each
 720 !      non-local projector: <p_lmn|c>
 721 !      ==============================================================
 722 
 723 !      Retrieve eventually <p_lmn|c> coeffs (and derivatives)
 724        if (cpopt>=2) then
 725          do ispinor=1,nspinor
 726            do ia=1,nincat
 727              gx(1:cplex,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)
 728            end do
 729          end do
 730        end if
 731        if (cpopt==4.and.ndgxdt>0) then
 732          ndgxdt_stored = cprjin(1,1)%ncpgr
 733          ishift=0
 734          if (((choice==2).or.(choice==3)).and.(ndgxdt_stored>ndgxdt).and.(signs==2)) ishift=idir-ndgxdt
 735          if ((choice==2).and.(ndgxdt_stored==9).and.(signs==2)) ishift=ishift+6
 736          if (choice==2.and.(ndgxdt_stored>ndgxdt).and.(signs==1)) ishift=ndgxdt_stored-ndgxdt
 737          if(cplex == 2) then
 738            do ispinor=1,nspinor
 739              do ia=1,nincat
 740                if (ndgxdt_stored==ndgxdt.or.(ndgxdt_stored>ndgxdt.and.((choice==2).or.(choice==3)))) then
 741                  dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)
 742                else if (signs==2.and.ndgxdt_stored==3) then
 743                  if (choice==5.or.choice==51.or.choice==52) then ! ndgxdt=1
 744                    dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir,1:nlmn)
 745                  else if (choice==8) then ! ndgxdt=2
 746                    idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
 747                    dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir1,1:nlmn)
 748                    dgxdt(1:2,2,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir2,1:nlmn)
 749                  else if (choice==81) then ! ndgxdt=1
 750                    idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
 751                    dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir2,1:nlmn)
 752                  end if
 753                end if
 754              end do
 755            end do
 756          else
 757            do ispinor=1,nspinor
 758              do ia=1,nincat
 759                do ilmn=1,nlmn
 760                  if (ndgxdt_stored==ndgxdt.or.(ndgxdt_stored>ndgxdt.and.((choice==2).or.(choice==3)))) then
 761                    do ii=1,ndgxdt
 762                      ic = cplex_dgxdt(ii)
 763                      dgxdt(1,ii,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)
 764                    end do
 765                  else if (signs==2.and.ndgxdt_stored==3) then
 766                    if (choice==5.or.choice==51.or.choice==52) then ! ndgxdt=1
 767                      dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir,ilmn)
 768                    else if (choice==8) then ! ndgxdt=2
 769                      idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
 770                      dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir1,ilmn)
 771                      dgxdt(1,2,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(2),idir2,ilmn)
 772                    else if (choice==81) then ! ndgxdt=1
 773                      idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
 774                      dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir2,ilmn)
 775                    end if
 776                  end if
 777                end do
 778              end do
 779            end do
 780          end if
 781        end if
 782 
 783 !      Computation or <p_lmn|c> (and derivatives) for this block of atoms
 784        if ((cpopt<4.and.choice_a/=-1).or.choice==8.or.choice==81) then
 785          call opernla_ylm(choice_a,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnlin,d2gxdt,dgxdt,ffnlin_typ,gx,&
 786 &         ia3,idir,indlmn_typ,istwf_k,kpgin_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpgin_,nlmn,&
 787 &         nloalg,npwin,nspinor,ph3din,signs,ucvol,vectin)
 788        end if
 789 
 790 !      Transfer result to output variable cprj (if requested)
 791 !      cprj(:)%cp receive the <p_i|Psi> factors (p_i: non-local projector)
 792 !      Be careful: cprj(:)%dcp does not exactly contain the derivative of cprj(:)%cp.
 793 !                  - Volume contributions (in case of strain derivative) are not included,
 794 !                  - Global coordinate transformation are not applied,
 795 !                  cprj(:)%dcp is meant to be a restart argument of the present nonlop routine.
 796        if (cpopt==0.or.cpopt==1) then
 797          do ispinor=1,nspinor
 798            do ia=1,nincat
 799              cprjin(iatm+ia,ispinor)%nlmn=nlmn
 800              cprjin(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)=gx(1:cplex,1:nlmn,ia,ispinor)
 801              if (cplex==1) cprjin(iatm+ia,ispinor)%cp(2,1:nlmn)=zero
 802            end do
 803          end do
 804        end if
 805        if ((cpopt==1.or.cpopt==3).and.ndgxdt>0) then
 806          ishift=0
 807          if ((choice==2).and.(cprjin(1,1)%ncpgr>ndgxdt)) ishift=cprjin(1,1)%ncpgr-ndgxdt
 808          if(cplex==2)then
 809            do ispinor=1,nspinor
 810              do ia=1,nincat
 811                cprjin(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)=dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)
 812 !               cprjin(iatm+ia,ispinor)%dcp(1:2,1:ndgxdt,1:nlmn)=dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)
 813              end do
 814            end do
 815          else
 816            do ispinor=1,nspinor
 817              do ia=1,nincat
 818                do ilmn=1,nlmn
 819                  do ii=1,ndgxdt
 820                    ic = cplex_dgxdt(ii) ; jc = 3 - ic
 821                    cprjin(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor)
 822                    cprjin(iatm+ia,ispinor)%dcp(jc,ii+ishift,ilmn)=zero
 823 !                   cprjin(iatm+ia,ispinor)%dcp(ic,ii,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor)
 824 !                   cprjin(iatm+ia,ispinor)%dcp(jc,ii,ilmn)=zero
 825                  end do
 826                end do
 827              end do
 828            end do
 829          end if
 830        end if
 831 
 832 !      If choice==0, that's all for these atoms !
 833        if (choice>0) then
 834          if(choice/=7) then
 835 !          Contraction from <p_i|c> to Sum_j[Dij.<p_j|c>] (and derivatives)
 836            call opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,dgxdt,dgxdtfac,dgxdtfac_sij,&
 837 &           d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,dimekbq,enl,gx,gxfac,gxfac_sij,&
 838 &           iatm,indlmn_typ,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,&
 839 &           nincat,nlmn,nspinor,nspinortot,optder,paw_opt,sij_typ)
 840          else
 841            gxfac_sij=gx
 842          end if
 843 
 844 !        Operate with the non-local potential on the projected scalars,
 845 !        in order to get contributions to energy/forces/stress/dyn.mat
 846 !        ==============================================================
 847          if (signs==1) then
 848            if (.not.present(cprjin_left)) then
 849              call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
 850 &             enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,nd2gxdt,ndgxdt,ndgxdtfac,&
 851 &             nincat,nlmn,nnlout,nspinor,paw_opt,strnlk)
 852            else
 853              ABI_ALLOCATE(gx_left,(cplex,nlmn,nincat,nspinor))
 854 !            Retrieve <p_lmn|c> coeffs
 855              do ispinor=1,nspinor
 856                do ia=1,nincat
 857                  gx_left(1:cplex,1:nlmn,ia,ispinor)=cprjin_left(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)
 858                end do
 859              end do
 860 !            TODO
 861 !            if (cpopt==4.and.ndgxdt>0) then
 862 !            do ispinor=1,nspinor
 863 !            do ia=1,nincat
 864 !            dgxdt_left(1:cplex,1:ndgxdt,1:nlmn,ia,ispinor)=cprjin_left(iatm+ia,ispinor)%dcp(1:cplex,1:ndgxdt,1:nlmn)
 865 !            end do
 866 !            end do
 867 !            end if
 868              call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
 869 &             enlk,enlout,fnlk,gx_left,gxfac,gxfac_sij,ia3,natom,nd2gxdt,ndgxdt,ndgxdtfac,&
 870 &             nincat,nlmn,nnlout,nspinor,paw_opt,strnlk)
 871              ABI_DEALLOCATE(gx_left)
 872            end if
 873          end if
 874 
 875 !        Operate with the non-local potential on the projected scalars,
 876 !        in order to get matrix element
 877 !        ==============================================================
 878          if (signs==2) then
 879 !          Prepare the phase factors if they were not already computed
 880            if(nloalg(2)<=0) then
 881              call ph1d3d(ia3,ia4,kgout,matblk,natom,npwout,n1,n2,n3,phkxredout,ph1d,ph3dout)
 882            end if
 883            call opernlb_ylm(choice_b,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,&
 884 &           d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnlout,ffnlout_typ,gxfac,gxfac_sij,ia3,&
 885 &           idir,indlmn_typ,kpgout_,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpgout_,nlmn,&
 886 &           nloalg,npwout,nspinor,paw_opt,ph3dout,svectout,ucvol,vectout)
 887          end if
 888 
 889        end if ! choice==0
 890 
 891 !      Deallocate temporary projected scalars
 892        ABI_DEALLOCATE(gx)
 893        ABI_DEALLOCATE(gxfac)
 894        ABI_DEALLOCATE(dgxdt)
 895        ABI_DEALLOCATE(dgxdtfac)
 896        ABI_DEALLOCATE(d2gxdt)
 897        ABI_DEALLOCATE(d2gxdtfac)
 898        ABI_DEALLOCATE(dgxdtfac_sij)
 899        ABI_DEALLOCATE(d2gxdtfac_sij)
 900        ABI_DEALLOCATE(gxfac_sij)
 901        ABI_DEALLOCATE(cplex_dgxdt)
 902        ABI_DEALLOCATE(cplex_d2gxdt)
 903 
 904 !      End sum on atom subset loop
 905        iatm=iatm+nincat;ia5=ia5+nincat
 906      end do
 907      !if (paw_opt>=2)  then
 908      ABI_DEALLOCATE(sij_typ)
 909      !end if
 910 
 911 !    End condition of existence of a non-local part
 912    else
 913      if (cpopt==0.or.cpopt==1) then
 914        do ispinor=1,nspinor
 915          do ia=1,nattyp(itypat)
 916            cprjin(iatm+ia,ispinor)%cp(:,1:nlmn)=zero
 917          end do
 918        end do
 919      end if
 920      if ((cpopt==1.or.cpopt==3).and.ndgxdt>0) then
 921        ishift=0
 922        if ((choice==2).and.(cprjin(1,1)%ncpgr>ndgxdt)) ishift=cprjin(1,1)%ncpgr-ndgxdt
 923        do ispinor=1,nspinor
 924          do ia=1,nattyp(itypat)
 925            cprjin(iatm+ia,ispinor)%dcp(:,1+ishift:ndgxdt+ishift,1:nlmn)=zero
 926 !           cprjin(iatm+ia,ispinor)%dcp(:,1:ndgxdt,1:nlmn)=zero
 927          end do
 928        end do
 929      end if
 930      iatm=iatm+nattyp(itypat)
 931    end if
 932 
 933 !  End atom type loop
 934    ia1=ia2+1
 935  end do
 936 
 937 
 938 !Reduction in case of parallelism
 939 !==============================================================
 940 
 941  if (signs==1.and.mpi_enreg%paral_spinor==1) then
 942    if (nnlout/=0) then
 943      call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr)
 944    end if
 945    if (choice==3.or.choice==6.or.choice==23) then
 946      call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr)
 947    end if
 948    if (choice==6) then
 949      call xmpi_sum(fnlk,mpi_enreg%comm_spinor,ierr)
 950      call xmpi_sum(strnlk,mpi_enreg%comm_spinor,ierr)
 951    end if
 952    if (choice==55) then
 953      call xmpi_sum(ddkk,mpi_enreg%comm_spinor,ierr)
 954    end if
 955  end if
 956 
 957 !Coordinate transformations
 958 !==============================================================
 959 
 960 !Need sometimes gmet
 961  if ((signs==1.and.paw_opt<=3).and. &
 962 & (choice==5 .or.choice==51.or.choice==52.or.choice==53.or.&
 963 & choice==54.or.choice==55)) then
 964    ABI_ALLOCATE(gmet,(3,3))
 965    gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
 966  end if
 967 
 968 !1st derivative wrt to strain (stress tensor):
 969 ! - convert from reduced to cartesian coordinates
 970 ! - substract volume contribution
 971  if ((choice==3.or.choice==23).and.signs==1.and.paw_opt<=3) then
 972    mu0=0 ! Shift to be applied in enlout array
 973    ABI_ALLOCATE(work1,(6))
 974    work1(1:6)=enlout(mu0+1:mu0+6)
 975    call strconv(work1,gprimd,work1)
 976    enlout(mu0+1:mu0+3)=(work1(1:3)-enlk)
 977    enlout(mu0+4:mu0+6)= work1(4:6)
 978    ABI_DEALLOCATE(work1)
 979  end if
 980 
 981 !1st derivative wrt to k wave vector (ddk):
 982 ! - convert from cartesian to reduced coordinates
 983  if ((choice==5.or.choice==53).and.signs==1.and.paw_opt<=3) then
 984    mu0=0 ! Shift to be applied in enlout array
 985    ABI_ALLOCATE(work1,(3))
 986    work1(:)=enlout(mu0+1:mu0+3)
 987    enlout(mu0+1:mu0+3)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3)
 988    ABI_DEALLOCATE(work1)
 989  end if
 990  if ((choice==51.or.choice==52).and.signs==1.and.paw_opt<=3) then
 991    mu0=0 ! Shift to be applied in enlout array
 992    ABI_ALLOCATE(work1,(3))
 993    do mu=1,2 ! Loop for Re,Im
 994      work1(1:3)=(/enlout(mu0+1),enlout(mu0+3),enlout(mu0+5)/)
 995      enlout(mu0+1)=gmet(1,1)*work1(1)+gmet(1,2)*work1(2)+gmet(1,3)*work1(3)
 996      enlout(mu0+3)=gmet(2,1)*work1(1)+gmet(2,2)*work1(2)+gmet(2,3)*work1(3)
 997      enlout(mu0+5)=gmet(3,1)*work1(1)+gmet(3,2)*work1(2)+gmet(3,3)*work1(3)
 998      mu0=mu0+1
 999    end do
1000    ABI_DEALLOCATE(work1)
1001  end if
1002 
1003 !2nd derivative wrt to k wave vector and atomic position (effective charges):
1004 ! - convert from cartesian to reduced coordinates
1005  if (choice==54.and.signs==1.and.paw_opt<=3) then
1006    mu0=0 ! Shift to be applied in enlout array
1007    ABI_ALLOCATE(work1,(3))
1008    ABI_ALLOCATE(work2,(3))
1009    do mu=1,3*natom
1010 !    First, real part
1011      work1(1)=enlout(mu0+1);work1(2)=enlout(mu0+3);work1(3)=enlout(mu0+5)
1012      work2(:)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3)
1013      enlout(mu0+1)=work2(1);enlout(mu0+3)=work2(2);enlout(mu0+5)=work2(3)
1014 !    Then imaginary part
1015      work1(1)=enlout(mu0+2);work1(2)=enlout(mu0+4);work1(3)=enlout(mu0+6)
1016      work2(:)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3)
1017      enlout(mu0+2)=work2(1);enlout(mu0+4)=work2(2);enlout(mu0+6)=work2(3)
1018      mu0=mu0+6
1019    end do
1020    ABI_DEALLOCATE(work1)
1021    ABI_DEALLOCATE(work2)
1022  end if
1023 
1024 !2nd derivative wrt to k wave vector and strain (piezoelectric tensor):
1025 ! - convert from cartesian to reduced coordinates (k point)
1026 ! - convert from reduced to cartesian coordinates (strain)
1027 ! - substract volume contribution
1028 ! - symetrize strain components
1029  if (choice==55.and.signs==1.and.paw_opt<=3) then
1030    ABI_ALLOCATE(work3,(2,3))
1031    ABI_ALLOCATE(work4,(2,3))
1032    ABI_ALLOCATE(work5,(2,3,6))
1033    ABI_ALLOCATE(work7,(2,3,6))
1034    ABI_ALLOCATE(work6,(2,3,3))
1035    do ic=1,3 ! gamma
1036      work5=zero
1037      do jc=1,3 ! nu
1038        do ii=1,3 ! lambda
1039          mu=(gamma(jc,ii)-1)*3+1
1040          work5(1,jc,ii)=gmet(ic,1)*enlout(2*mu-1)+gmet(ic,2)*enlout(2*mu+1) &
1041 &         +gmet(ic,3)*enlout(2*mu+3)
1042          work5(2,jc,ii)=gmet(ic,1)*enlout(2*mu  )+gmet(ic,2)*enlout(2*mu+2) &
1043 &         +gmet(ic,3)*enlout(2*mu+4)
1044        end do
1045      end do
1046      work6=zero
1047      do jc=1,3 ! nu
1048        do ii=1,3 ! beta
1049          work6(1:cplex,ii,jc)=gprimd(ii,1)*work5(1:cplex,jc,1)+gprimd(ii,2)*work5(1:cplex,jc,2) &
1050 &         +gprimd(ii,3)*work5(1:cplex,jc,3)
1051        end do
1052      end do
1053      do jc=1,3 ! alpha
1054        do ii=1,3 ! beta
1055          mu=gamma(jc,ii)
1056          work7(1:cplex,ic,mu)=gprimd(jc,1)*work6(1:cplex,ii,1)+gprimd(jc,2)*work6(1:cplex,ii,2) &
1057 &         +gprimd(jc,3)*work6(1:cplex,ii,3)
1058        end do
1059      end do
1060    end do ! gamma
1061 
1062    do ii=1,3 ! alpha
1063      work3(1,ii)=gprimd(ii,1)*ddkk(2*1-1)+gprimd(ii,2)*ddkk(2*2-1) &
1064 &     +gprimd(ii,3)*ddkk(2*3-1)
1065      work3(2,ii)=gprimd(ii,1)*ddkk(2*1  )+gprimd(ii,2)*ddkk(2*2  ) &
1066 &     +gprimd(ii,3)*ddkk(2*3  )
1067    end do
1068    do ii=1,3 ! gamma
1069      work4(1,ii)=gmet(ii,1)*ddkk(2*1-1)+gmet(ii,2)*ddkk(2*2-1) &
1070 &     +gmet(ii,3)*ddkk(2*3-1)
1071      work4(2,ii)=gmet(ii,1)*ddkk(2*1  )+gmet(ii,2)*ddkk(2*2  ) &
1072 &     +gmet(ii,3)*ddkk(2*3  )
1073    end do
1074 
1075    do mu=1,6
1076      ii=alpha(mu) ! alpha
1077      ic=beta(mu) ! beta
1078      do jc=1,3 ! gamma
1079        work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-half &
1080 &       *(gprimd(ic,jc)*work3(1:cplex,ii)+gprimd(ii,jc)*work3(1:cplex,ic))
1081        if (ii==ic) work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-work4(1:cplex,jc)
1082      end do
1083    end do
1084    do mu=1,6 ! alpha,beta
1085      do nu=1,3 ! gamma
1086        mu0=3*(mu-1)+nu
1087        enlout(2*mu0-1)=work7(1,nu,mu)
1088        enlout(2*mu0  )=work7(2,nu,mu)
1089      end do
1090    end do
1091    ABI_DEALLOCATE(gmet)
1092    ABI_DEALLOCATE(work3)
1093    ABI_DEALLOCATE(work4)
1094    ABI_DEALLOCATE(work5)
1095    ABI_DEALLOCATE(work6)
1096    ABI_DEALLOCATE(work7)
1097  end if
1098 
1099 !2nd derivative wrt to 2 k wave vectors (effective mass):
1100 ! - convert from cartesian to reduced coordinates
1101  if ((choice==8.or.choice==81).and.signs==1.and.paw_opt<=3) then
1102    mu0=0 ! Shift to be applied in enlout array
1103    ABI_ALLOCATE(work3,(3,3))
1104    ABI_ALLOCATE(work4,(3,3))
1105    mua=1;if (choice==81) mua=2
1106    do ii=1,mua ! Loop Re,Im
1107      if (choice==8) then ! enlout is real in Voigt notation
1108        work3(1,1)=enlout(mu0+1) ; work3(1,2)=enlout(mu0+6) ; work3(1,3)=enlout(mu0+5)
1109        work3(2,1)=enlout(mu0+6) ; work3(2,2)=enlout(mu0+2) ; work3(2,3)=enlout(mu0+4)
1110        work3(3,1)=enlout(mu0+5) ; work3(3,2)=enlout(mu0+4) ; work3(3,3)=enlout(mu0+3)
1111      else                ! enlout is complex in matrix notation
1112        work3(1,1)=enlout(mu0+1 ) ; work3(1,2)=enlout(mu0+3 ) ; work3(1,3)=enlout(mu0+5 )
1113        work3(2,1)=enlout(mu0+7 ) ; work3(2,2)=enlout(mu0+9 ) ; work3(2,3)=enlout(mu0+11)
1114        work3(3,1)=enlout(mu0+13) ; work3(3,2)=enlout(mu0+15) ; work3(3,3)=enlout(mu0+17)
1115      end if
1116      do mu=1,3
1117        work4(:,mu)=gprimd(:,1)*work3(mu,1)+gprimd(:,2)*work3(mu,2)+gprimd(:,3)*work3(mu,3)
1118      end do
1119      do mu=1,3
1120        work3(:,mu)=gprimd(:,1)*work4(mu,1)+gprimd(:,2)*work4(mu,2)+gprimd(:,3)*work4(mu,3)
1121      end do
1122      do mu=1,3
1123        work4(:,mu)=gprimd(1,:)*work3(mu,1)+gprimd(2,:)*work3(mu,2)+gprimd(3,:)*work3(mu,3)
1124      end do
1125      do mu=1,3
1126        work3(:,mu)=gprimd(1,:)*work4(mu,1)+gprimd(2,:)*work4(mu,2)+gprimd(3,:)*work4(mu,3)
1127      end do
1128      if (choice==8) then ! enlout is real in Voigt notation
1129        enlout(mu0+1) = work3(1,1) ; enlout(mu0+2) = work3(2,2) ; enlout(mu0+3) = work3(3,3)
1130        enlout(mu0+4) = work3(3,2) ; enlout(mu0+5) = work3(1,3) ; enlout(mu0+6) = work3(2,1)
1131      else                ! enlout is complex in matrix notation
1132        enlout(mu0+1 )=work3(1,1) ; enlout(mu0+3 )=work3(1,2) ; enlout(mu0+5 )=work3(1,3)
1133        enlout(mu0+7 )=work3(2,1) ; enlout(mu0+9 )=work3(2,2) ; enlout(mu0+11)=work3(2,3)
1134        enlout(mu0+13)=work3(3,1) ; enlout(mu0+15)=work3(3,2) ; enlout(mu0+17)=work3(3,3)
1135      end if
1136      mu0=mu0+1
1137    end do
1138    ABI_DEALLOCATE(work3)
1139    ABI_DEALLOCATE(work4)
1140  end if
1141 
1142 !2nd derivative wrt to 2 strains (elastic tensor):
1143 ! - convert from reduced to cartesian coordinates
1144 ! - substract volume contribution
1145  if (choice==6.and.signs==1.and.paw_opt<=3) then
1146    mu0=0 ! Shift to be applied in enlout array
1147    ABI_ALLOCATE(work1,(6))
1148    ABI_ALLOCATE(work2,(6))
1149    ABI_ALLOCATE(work3,(6+3*natom,6))
1150    work3(:,:)=reshape(enlout(mu0+1:mu0+6*(6+3*natom)),(/6+3*natom,6/))
1151    do mu=1,6
1152      call strconv(work3(1:6,mu),gprimd,work3(1:6,mu))
1153    end do
1154    do mu=1,6+3*natom
1155      work1(1:6)=work3(mu,1:6)
1156      call strconv(work1,gprimd,work2)
1157      work3(mu,1:6)=work2(1:6)
1158    end do
1159    enlout(mu0+1:mu0+6*(6+3*natom))=reshape(work3(:,:),(/6*(6+3*natom)/))
1160    ABI_DEALLOCATE(work1)
1161    ABI_DEALLOCATE(work2)
1162    ABI_DEALLOCATE(work3)
1163    call strconv(strnlk,gprimd,strnlk)
1164    do mub=1,6
1165      nub1=alpha(mub);nub2=beta(mub)
1166      do mua=1,6
1167        mu=mu0+mua+(3*natom+6)*(mub-1)
1168        nua1=alpha(mua);nua2=beta(mua)
1169        if (mua<=3.and.mub<=3) enlout(mu)=enlout(mu)+enlk
1170        if (mua<=3) enlout(mu)=enlout(mu)-strnlk(mub)
1171        if (mub<=3) enlout(mu)=enlout(mu)-strnlk(mua)
1172        if (nub1==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub2))
1173        if (nub2==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub1))
1174        if (nub1==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub2))
1175        if (nub2==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub1))
1176      end do
1177      if (mub<=3) then
1178        do nua1=1,natom
1179          nua2=3*(nua1-1);mu=mu0+nua2+6+(3*natom+6)*(mub-1)
1180          enlout(mu+1:mu+3)=enlout(mu+1:mu+3)-fnlk(nua2+1:nua2+3)
1181        end do
1182      end if
1183    end do
1184  end if
1185 
1186  if (allocated(gmet)) then
1187    ABI_DEALLOCATE(gmet)
1188  end if
1189 
1190 !Final deallocations
1191 !==============================================================
1192 
1193  if (signs==1)  then
1194    ABI_DEALLOCATE(fnlk)
1195    ABI_DEALLOCATE(ddkk)
1196    ABI_DEALLOCATE(strnlk)
1197  end if
1198  if (nkpgin<nkpgin_) then
1199    ABI_DEALLOCATE(kpgin_)
1200  end if
1201  if (nkpgout<nkpgout_) then
1202    ABI_DEALLOCATE(kpgout_)
1203  end if
1204 
1205  DBG_EXIT("COLL")
1206 
1207 end subroutine nonlop_ylm