TABLE OF CONTENTS


ABINIT/m_opernlc_ylm [ Modules ]

[ Top ] [ Modules ]

NAME

  m_opernlc_ylm

FUNCTION

COPYRIGHT

  Copyright (C) 2008-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_opernlc_ylm
27 
28  use defs_basis
29  use defs_abitypes
30  use m_errors
31  use m_abicore
32  use m_xmpi
33 
34  implicit none
35 
36  private

ABINIT/opernlc_ylm [ Functions ]

[ Top ] [ Functions ]

NAME

 opernlc_ylm

FUNCTION

 * Operate with the non-local part of the hamiltonian,
   in order to reduce projected scalars
 * Operate with the non-local projectors and the overlap matrix,
   in order to reduce projected scalars

INPUTS

  atindx1(natom)=index table for atoms (gives the absolute index of
                 an atom from its rank in a block of atoms)
  cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
        2 if <p_lmn|c> scalars are complex
  cplex_dgxdt(ndgxdt) = used only when cplex = 1
             cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:)   is real, 2 if it is pure imaginary
  cplex_enl=1 if enl factors are real, 2 if they are complex
  cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
  dgxdt(cplex,ndgxdt,nlmn,nincat)=grads of projected scalars (only if optder>0)
  dimenl1,dimenl2=dimensions of enl (see enl)
  dimekbq=1 if enl factors do not contain a exp(-iqR) phase, 2 is they do
  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
  gx(cplex,nlmn,nincat*abs(enl_opt))= projected scalars
  iatm=absolute rank of first atom of the current block of atoms
  indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
  itypat=type of atoms
  lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
  mpi_enreg=informations about MPI parallelization
  natom=number of atoms in cell
  ndgxdt=second dimension of dgxdt
  ndgxdtfac=second dimension of dgxdtfac
  nincat=number of atoms in the subset here treated
  nlmn=number of (l,m,n) numbers for current type of atom
  nspinor= number of spinorial components of the wavefunctions (on current proc)
  nspinortot=total number of spinorial components of the wavefunctions
  optder=0=only gxfac is computed, 1=both gxfac and dgxdtfac are computed
         2=gxfac, dgxdtfac and d2gxdtfac are computed
  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)
  sij(nlm*(nlmn+1)/2)=overlap matrix components (only if paw_opt=2, 3 or 4)

OUTPUT

  if (paw_opt=0, 1, 2 or 4)
    gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
  if (paw_opt=3 or 4)
    gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
  if (optder==1.and.paw_opt=0, 1, 2 or 4)
    dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Vnl (NL operator)
  if (optder==1.and.paw_opt=3 or 4)
    dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Sij (overlap)

NOTES

 This routine operates for one type of atom, and within this given type of atom,
 for a subset of at most nincat atoms.

 About the non-local factors symmetry:
   - The lower triangular part of the Dij matrix can be deduced from the upper one
     with the following relation: D^s2s1_ji = (D^s1s2_ij)^*
     where s1,s2 are spinor components
   - The Dij factors can contain a exp(-iqR) phase
     This phase does not have to be included in the symmetry rule
     For that reason, we first apply the real part (cos(qR).D^s1s2_ij)
     then, we apply the imaginary part (-sin(qR).D^s1s2_ij)

PARENTS

      m_gemm_nonlop,nonlop_ylm

CHILDREN

      xmpi_sum

SOURCE

 135 subroutine opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,&
 136 &          dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,dimekbq,enl,&
 137 &          gx,gxfac,gxfac_sij,iatm,indlmn,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,&
 138 &          nd2gxdt,nd2gxdtfac,nincat,nlmn,nspinor,nspinortot,optder,paw_opt,sij)
 139 
 140 
 141 !This section has been created automatically by the script Abilint (TD).
 142 !Do not modify the following lines by hand.
 143 #undef ABI_FUNC
 144 #define ABI_FUNC 'opernlc_ylm'
 145 !End of the abilint section
 146 
 147  implicit none
 148 
 149 !Arguments ------------------------------------
 150 !scalars
 151  integer,intent(in) :: cplex,cplex_enl,cplex_fac,dimenl1,dimenl2,dimekbq,iatm,itypat
 152  integer,intent(in) :: natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,nincat,nspinor,nspinortot,optder,paw_opt
 153  integer,intent(inout) :: nlmn
 154  real(dp) :: lambda
 155  type(MPI_type) , intent(in) :: mpi_enreg
 156 !arrays
 157  integer,intent(in) :: atindx1(natom),indlmn(6,nlmn),cplex_dgxdt(ndgxdt),cplex_d2gxdt(nd2gxdt)
 158  real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)
 159  real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)
 160  real(dp),intent(in),target :: enl(dimenl1,dimenl2,nspinortot**2,dimekbq)
 161  real(dp),intent(inout) :: gx(cplex,nlmn,nincat,nspinor)
 162  real(dp),intent(in) :: sij(((paw_opt+1)/3)*nlmn*(nlmn+1)/2)
 163  real(dp),intent(out),target :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)
 164  real(dp),intent(out) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3))
 165  real(dp),intent(out),target :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor)
 166  real(dp),intent(out) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat,nspinor*(paw_opt/3))
 167  real(dp),intent(out),target :: gxfac(cplex_fac,nlmn,nincat,nspinor)
 168  real(dp),intent(out) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
 169 
 170 !Local variables-------------------------------
 171 !Arrays
 172 !scalars
 173  integer :: cplex_,ia,ierr,ijlmn,ijspin,ilm,ilmn,i0lmn,iln,index_enl,iphase,ispinor,ispinor_index
 174  integer :: j0lmn,jilmn,jispin,jjlmn,jlm,jlmn,jspinor,jspinor_index,mu,shift
 175  real(dp) :: sijr
 176 !arrays
 177  real(dp) :: enl_(2),gxfi(2),gxi(cplex),gxj(cplex)
 178  real(dp),allocatable :: d2gxdtfac_offdiag(:,:,:,:,:),dgxdtfac_offdiag(:,:,:,:,:)
 179  real(dp),allocatable :: gxfac_offdiag(:,:,:,:),gxfj(:,:)
 180  real(dp),pointer :: d2gxdtfac_(:,:,:,:,:),dgxdtfac_(:,:,:,:,:),gxfac_(:,:,:,:)
 181  real(dp),pointer :: enl_ptr(:,:,:)
 182 
 183 ! *************************************************************************
 184 
 185  DBG_ENTER("COLL")
 186 
 187 !Parallelization over spinors treatment
 188  shift=0;if (mpi_enreg%paral_spinor==1) shift=mpi_enreg%me_spinor
 189 
 190 !When Enl factors contain a exp(-iqR) phase:
 191 ! - We loop over the real and imaginary parts
 192 ! - We need an additional memory space
 193  do iphase=1,dimekbq
 194   if (paw_opt==3) cycle
 195   if (iphase==1) then
 196    gxfac_ => gxfac ; dgxdtfac_ => dgxdtfac ; d2gxdtfac_ => d2gxdtfac
 197   else
 198    ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac==1 when dimekbq=2!")
 199    ABI_ALLOCATE(gxfac_,(cplex_fac,nlmn,nincat,nspinor))
 200    if (optder>=1) then
 201     ABI_ALLOCATE(dgxdtfac_,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor))
 202    end if
 203    if (optder>=2) then
 204     ABI_ALLOCATE(d2gxdtfac_,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor))
 205    end if
 206   end if
 207   gxfac_=zero
 208   if (optder>=1) dgxdtfac_=zero
 209   if (optder>=2) d2gxdtfac_=zero
 210   enl_ptr => enl(:,:,:,iphase)
 211 
 212   
 213 !Accumulate gxfac related to non-local operator (Norm-conserving)
 214 !-------------------------------------------------------------------
 215   if (paw_opt==0) then
 216    !Enl is E(Kleinman-Bylander)
 217    ABI_CHECK(cplex_enl/=2,"BUG: invalid cplex_enl=2!")
 218    ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!")
 219 !$OMP PARALLEL &
 220 !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_)
 221 !$OMP DO COLLAPSE(3)
 222    do ispinor=1,nspinor
 223      do ia=1,nincat
 224        do ilmn=1,nlmn
 225          ispinor_index=ispinor+shift
 226          iln=indlmn(5,ilmn)
 227          enl_(1)=enl_ptr(iln,itypat,ispinor_index)
 228          gxfac_(1:cplex,ilmn,ia,ispinor)=enl_(1)*gx(1:cplex,ilmn,ia,ispinor)
 229        end do
 230      end do
 231    end do
 232 !$OMP END DO
 233 !$OMP END PARALLEL
 234   end if
 235 
 236 !Accumulate gxfac related to nonlocal operator (PAW)
 237 !-------------------------------------------------------------------
 238   if (paw_opt==1.or.paw_opt==2.or.paw_opt==4) then
 239    !Enl is psp strength Dij or (Dij-lambda.Sij)
 240 
 241 !  === Diagonal term(s) (up-up, down-down)
 242 
 243 !  1-Enl is real
 244    if (cplex_enl==1) then
 245 !$OMP PARALLEL &
 246 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
 247 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,ijlmn,gxi)
 248 !$OMP DO COLLAPSE(3)
 249      do ispinor=1,nspinor
 250        do ia=1,nincat
 251          do jlmn=1,nlmn
 252            ispinor_index=ispinor+shift
 253            index_enl=atindx1(iatm+ia)
 254            j0lmn=jlmn*(jlmn-1)/2
 255            jjlmn=j0lmn+jlmn
 256            enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
 257            if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 258            gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
 259            gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxj(1:cplex)
 260            do ilmn=1,jlmn-1
 261              ijlmn=j0lmn+ilmn
 262              enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
 263              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 264              gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
 265              gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex)
 266 #if !defined HAVE_OPENMP
 267              gxfac_(1:cplex,ilmn,ia,ispinor)=gxfac_(1:cplex,ilmn,ia,ispinor)+enl_(1)*gxj(1:cplex)
 268 #endif
 269            end do
 270 #if defined HAVE_OPENMP
 271            if(jlmn<nlmn) then
 272              do ilmn=jlmn+1,nlmn
 273                i0lmn=(ilmn*(ilmn-1)/2)
 274                ijlmn=i0lmn+jlmn
 275                enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
 276                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 277                gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
 278                gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex)
 279              end do
 280            end if
 281 #endif
 282          end do
 283        end do
 284      end do
 285 !$OMP END DO
 286 !$OMP END PARALLEL
 287 
 288 !    2-Enl is complex  ===== D^ss'_ij=D^s's_ji^*
 289    else
 290      ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
 291 
 292      if (nspinortot==1) then ! -------------> NO SPINORS
 293 
 294 !$OMP PARALLEL &
 295 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,ijlmn,gxi)
 296        do ia=1,nincat
 297          index_enl=atindx1(iatm+ia)
 298 !$OMP DO
 299          do jlmn=1,nlmn
 300            j0lmn=jlmn*(jlmn-1)/2
 301            jjlmn=j0lmn+jlmn
 302            enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
 303            if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 304            gxj(1:cplex)=gx(1:cplex,jlmn,ia,1)
 305            gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxj(1)
 306            if (cplex==2) gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxj(2)
 307            do ilmn=1,jlmn-1
 308              ijlmn=j0lmn+ilmn
 309              enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
 310              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 311              gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
 312              gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxi(1)
 313              gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)-enl_(2)*gxi(1)
 314 #if !defined HAVE_OPENMP
 315              gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1)+enl_(1)*gxj(1)
 316              gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1)+enl_(2)*gxj(1)
 317 #endif
 318              if (cplex==2) then
 319                gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(2)*gxi(2)
 320                gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxi(2)
 321 #if !defined HAVE_OPENMP
 322                gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1)-enl_(2)*gxj(2)
 323                gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1)+enl_(1)*gxj(2)
 324 #endif
 325              end if
 326            end do
 327 #if defined HAVE_OPENMP
 328            if(jlmn<nlmn) then
 329              do ilmn=jlmn+1,nlmn
 330                i0lmn=ilmn*(ilmn-1)/2
 331                ijlmn=i0lmn+jlmn
 332                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
 333                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 334                gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
 335                gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxi(1)
 336                gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(2)*gxi(1)
 337                if (cplex==2) then
 338                  gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)-enl_(2)*gxi(2)
 339                  gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxi(2)
 340                end if
 341              end do
 342            end if
 343 #endif
 344          end do
 345 !$OMP END DO
 346        end do
 347 !$OMP END PARALLEL
 348 
 349      else ! -------------> SPINORIAL CASE
 350 
 351 !  === Diagonal term(s) (up-up, down-down)
 352 
 353 !$OMP PARALLEL &
 354 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
 355 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,ijlmn,gxi)
 356        do ispinor=1,nspinor
 357          ispinor_index=ispinor+shift
 358          do ia=1,nincat
 359            index_enl=atindx1(iatm+ia)
 360 !$OMP DO
 361            do jlmn=1,nlmn
 362              j0lmn=jlmn*(jlmn-1)/2
 363              jjlmn=j0lmn+jlmn
 364              enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
 365              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 366              gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
 367              gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxj(1)
 368              if (cplex==2)  gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxj(2)
 369              do ilmn=1,jlmn-1
 370                ijlmn=j0lmn+ilmn
 371                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
 372                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 373                gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
 374                gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
 375                gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)-enl_(2)*gxi(1)
 376 #if !defined HAVE_OPENMP
 377                gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)+enl_(1)*gxj(1)
 378                gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(2)*gxj(1)
 379 #endif
 380                if (cplex==2) then
 381                  gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(2)*gxi(2)
 382                  gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
 383 #if !defined HAVE_OPENMP
 384                  gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)-enl_(2)*gxj(2)
 385                  gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(1)*gxj(2)
 386 #endif
 387                end if
 388              end do
 389 #if defined HAVE_OPENMP
 390              if(jlmn<nlmn) then
 391                do ilmn=jlmn+1,nlmn
 392                  i0lmn=ilmn*(ilmn-1)/2
 393                  ijlmn=i0lmn+jlmn
 394                  enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
 395                  if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 396                  gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
 397                  gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
 398                  gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(2)*gxi(1)
 399                  if (cplex==2) then
 400                    gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)-enl_(2)*gxi(2)
 401                    gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
 402                  end if
 403                end do
 404              end if
 405 #endif
 406            end do
 407 !$OMP END DO
 408          end do
 409        end do
 410 !$OMP END PARALLEL
 411      end if !nspinortot
 412    end if !complex_enl
 413 
 414 !  === Off-diagonal term(s) (up-down, down-up)
 415 
 416 !  --- No parallelization over spinors ---
 417    if (nspinortot==2.and.nspinor==nspinortot) then
 418      ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
 419      ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex)!")
 420 !$OMP PARALLEL &
 421 !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), &
 422 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxi,gxj,ilmn,ijlmn)
 423      do ispinor=1,nspinortot
 424        jspinor=3-ispinor
 425        do ia=1,nincat
 426          index_enl=atindx1(iatm+ia)
 427 !$OMP DO
 428          do jlmn=1,nlmn
 429            j0lmn=jlmn*(jlmn-1)/2
 430            jjlmn=j0lmn+jlmn
 431            enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor )
 432            gxi(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
 433            gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(1)*gxi(1)
 434            gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)-enl_(2)*gxi(1)
 435            if (cplex==2) then
 436              gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(2)*gxi(2)
 437              gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)+enl_(1)*gxi(2)
 438            end if
 439 #if !defined HAVE_OPENMP
 440            gxj(1:cplex)=gx(1:cplex,jlmn,ia,jspinor)
 441 #endif
 442            do ilmn=1,jlmn-1
 443              ijlmn=j0lmn+ilmn
 444              enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
 445              gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
 446              gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(1)*gxi(1)
 447              gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)-enl_(2)*gxi(1)
 448 #if !defined HAVE_OPENMP
 449              gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)+enl_(1)*gxj(1)
 450              gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(2)*gxj(1)
 451 #endif
 452              if (cplex==2) then
 453                gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(2)*gxi(2)
 454                gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)+enl_(1)*gxi(2)
 455 #if !defined HAVE_OPENMP
 456                gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)-enl_(2)*gxj(2)
 457                gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(1)*gxj(2)
 458 #endif
 459              end if
 460            end do
 461 #if defined HAVE_OPENMP
 462            if(jlmn<nlmn) then
 463              do ilmn=jlmn+1,nlmn
 464                i0lmn=ilmn*(ilmn-1)/2
 465                ijlmn=i0lmn+jlmn
 466                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
 467                gxi(1:cplex)=gx(1:cplex,ilmn,ia,jspinor)
 468                gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
 469                gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(2)*gxi(1)
 470                if (cplex==2) then
 471                  gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)-enl_(2)*gxi(2)
 472                  gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
 473                end if
 474              end do
 475            end if
 476 #endif
 477          end do
 478 !$OMP END DO
 479        end do
 480      end do
 481 !$OMP END PARALLEL
 482 
 483 !    --- Parallelization over spinors ---
 484    else if (nspinortot==2.and.nspinor/=nspinortot) then
 485      ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
 486      ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
 487      ABI_ALLOCATE(gxfac_offdiag,(cplex_fac,nlmn,nincat,nspinortot))
 488 !$OMP WORKSHARE
 489      gxfac_offdiag(:,:,:,:)=zero
 490 !$OMP END WORKSHARE
 491      ispinor_index=mpi_enreg%me_spinor+1
 492      jspinor_index=3-ispinor_index
 493      if (ispinor_index==1) then
 494        ijspin=3;jispin=4
 495      else
 496        ijspin=4;jispin=3
 497      end if
 498 !$OMP PARALLEL &
 499 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,gxi)
 500      do ia=1,nincat
 501        index_enl=atindx1(iatm+ia)
 502 !$OMP DO
 503        do jlmn=1,nlmn
 504          j0lmn=jlmn*(jlmn-1)/2
 505          do ilmn=1,nlmn
 506            i0lmn=ilmn*(ilmn-1)/2
 507            if (ilmn<=jlmn) then
 508              ijlmn=j0lmn+ilmn
 509              enl_(1)= enl_ptr(2*ijlmn-1,index_enl,ijspin)
 510              enl_(2)=-enl_ptr(2*ijlmn  ,index_enl,ijspin)
 511            else
 512              jilmn=i0lmn+jlmn
 513              enl_(1)= enl_ptr(2*jilmn-1,index_enl,jispin)
 514              enl_(2)= enl_ptr(2*jilmn  ,index_enl,jispin)
 515            end if
 516            gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
 517            gxfac_offdiag(1,jlmn,ia,jspinor_index)= &
 518 &            gxfac_offdiag(1,jlmn,ia,jspinor_index)+enl_(1)*gxi(1)
 519            gxfac_offdiag(2,jlmn,ia,jspinor_index)= &
 520 &            gxfac_offdiag(2,jlmn,ia,jspinor_index)+enl_(2)*gxi(1)
 521            if (cplex==2) then
 522              gxfac_offdiag(1,jlmn,ia,jspinor_index)= &
 523 &              gxfac_offdiag(1,jlmn,ia,jspinor_index)-enl_(2)*gxi(2)
 524              gxfac_offdiag(2,jlmn,ia,jspinor_index)= &
 525 &              gxfac_offdiag(2,jlmn,ia,jspinor_index)+enl_(1)*gxi(2)
 526            end if
 527          end do !ilmn
 528        end do !jlmn
 529 !$OMP END DO
 530      end do !iat
 531 !$OMP END PARALLEL
 532      call xmpi_sum(gxfac_offdiag,mpi_enreg%comm_spinor,ierr)
 533      gxfac_(:,:,:,1)=gxfac_(:,:,:,1)+gxfac_offdiag(:,:,:,ispinor_index)
 534      ABI_DEALLOCATE(gxfac_offdiag)
 535    end if
 536 
 537   end if !paw_opt
 538 
 539 !Accumulate dgxdtfac related to nonlocal operator (Norm-conserving)
 540 !-------------------------------------------------------------------
 541   if (optder>=1.and.paw_opt==0) then
 542    !Enl is E(Kleinman-Bylander)
 543    ABI_CHECK(cplex_enl==1,"BUG: invalid cplex_enl/=1!")
 544    ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!")
 545 !$OMP PARALLEL &
 546 !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_,mu)
 547    do ispinor=1,nspinor
 548      ispinor_index = ispinor + shift
 549      do ia=1,nincat
 550 !$OMP DO
 551        do ilmn=1,nlmn
 552          iln=indlmn(5,ilmn)
 553          enl_(1)=enl_ptr(iln,itypat,ispinor_index)
 554          do mu=1,ndgxdtfac
 555            dgxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=enl_(1)*dgxdt(1:cplex,mu,ilmn,ia,ispinor)
 556          end do
 557        end do
 558 !$OMP END DO
 559      end do
 560    end do
 561 !$OMP END PARALLEL
 562   end if
 563 
 564 !Accumulate dgxdtfac related to nonlocal operator (PAW)
 565 !-------------------------------------------------------------------
 566   if (optder>=1.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then
 567    !Enl is psp strength Dij or (Dij-lambda.Sij)
 568 
 569 !  === Diagonal term(s) (up-up, down-down)
 570 
 571 !  1-Enl is real
 572    if (cplex_enl==1) then
 573 !$OMP PARALLEL &
 574 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi)
 575      ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac))
 576      do ispinor=1,nspinor
 577        ispinor_index=ispinor+shift
 578        do ia=1,nincat
 579          index_enl=atindx1(iatm+ia)
 580 !$OMP DO
 581          do jlmn=1,nlmn
 582            j0lmn=jlmn*(jlmn-1)/2
 583            jjlmn=j0lmn+jlmn
 584            enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
 585            if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 586            do mu=1,ndgxdtfac
 587              gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
 588              dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
 589            end do
 590            do ilmn=1,jlmn-1
 591              ijlmn=j0lmn+ilmn
 592              enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
 593              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 594              do mu=1,ndgxdtfac
 595                gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
 596                dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
 597 #if !defined HAVE_OPENMP
 598                dgxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
 599 #endif
 600              end do
 601            end do
 602 #if defined HAVE_OPENMP
 603            if(jlmn<nlmn) then
 604              do ilmn=jlmn+1,nlmn
 605                i0lmn=ilmn*(ilmn-1)/2
 606                ijlmn=i0lmn+jlmn
 607                enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
 608                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 609                do mu=1,ndgxdtfac
 610                  gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
 611                  dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
 612                end do
 613              end do
 614            end if
 615 #endif
 616          end do
 617 !$OMP END DO
 618        end do
 619      end do
 620      ABI_DEALLOCATE(gxfj)
 621 !$OMP END PARALLEL
 622 
 623 !    2-Enl is complex  ===== D^ss'_ij=D^s's_ji^*
 624    else
 625      ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
 626 
 627      if (nspinortot==1) then ! -------------> NO SPINORS
 628 
 629 !$OMP PARALLEL &
 630 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi)
 631        ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac))
 632        do ia=1,nincat
 633          index_enl=atindx1(iatm+ia)
 634 !$OMP DO
 635          do jlmn=1,nlmn
 636            j0lmn=jlmn*(jlmn-1)/2
 637            jjlmn=j0lmn+jlmn
 638            enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
 639            if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 640            do mu=1,ndgxdtfac
 641              if(cplex_dgxdt(mu)==2)then
 642                cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,1)
 643              else
 644                cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,1)
 645              end if
 646              dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu)
 647              if (cplex_==2) dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu)
 648            end do
 649            do ilmn=1,jlmn-1
 650              ijlmn=j0lmn+ilmn
 651              enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
 652              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 653              do mu=1,ndgxdtfac
 654                if(cplex_dgxdt(mu)==2)then
 655                  cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
 656                else
 657                  cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
 658                end if
 659                dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
 660                dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)-enl_(2)*gxfi(1)
 661 #if !defined HAVE_OPENMP
 662                dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu)
 663                dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu)
 664 #endif
 665                if (cplex_==2) then
 666                  dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(2)*gxfi(2)
 667                  dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
 668 #if !defined HAVE_OPENMP
 669                  dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu)
 670                  dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1)+enl_(1)*gxfj(2,mu)
 671 #endif
 672                end if
 673              end do
 674            end do
 675 #if defined HAVE_OPENMP
 676            if(jlmn<nlmn) then
 677              do ilmn=jlmn+1,nlmn
 678                i0lmn=ilmn*(ilmn-1)/2
 679                ijlmn=i0lmn+jlmn
 680                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
 681                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 682                do mu=1,ndgxdtfac
 683                  if(cplex_dgxdt(mu)==2)then
 684                    cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
 685                  else
 686                    cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
 687                  end if
 688                  dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
 689                  dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1)
 690                  if (cplex_==2) then
 691                    dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2)
 692                    dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
 693                  end if
 694                end do
 695              end do
 696            end if
 697 #endif
 698          end do
 699 !$OMP END DO
 700        end do
 701        ABI_DEALLOCATE(gxfj)
 702 !$OMP END PARALLEL
 703 
 704      else ! -------------> SPINORIAL CASE
 705 
 706 !  === Diagonal term(s) (up-up, down-down)
 707 
 708 !$OMP PARALLEL &
 709 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
 710 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi)
 711        ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac))
 712        do ispinor=1,nspinor
 713          ispinor_index = ispinor + shift
 714          do ia=1,nincat
 715            index_enl=atindx1(iatm+ia)
 716 !$OMP DO
 717            do jlmn=1,nlmn
 718              j0lmn=jlmn*(jlmn-1)/2
 719              jjlmn=j0lmn+jlmn
 720              enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
 721              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 722              do mu=1,ndgxdtfac
 723                if(cplex_dgxdt(mu)==2)then
 724                  cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,ispinor)
 725                else
 726                  cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
 727                end if
 728                dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
 729                if (cplex_==2) dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
 730              end do
 731              do ilmn=1,jlmn-1
 732                ijlmn=j0lmn+ilmn
 733                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
 734                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 735                do mu=1,ndgxdtfac
 736                  if(cplex_dgxdt(mu)==2)then
 737                    cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor)
 738                  else
 739                    cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
 740                  end if
 741                  dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
 742                  dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(1)
 743 #if !defined HAVE_OPENMP
 744                  dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
 745                  dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
 746 #endif
 747                  if (cplex_==2) then
 748                    dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(2)
 749                    dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
 750 #if !defined HAVE_OPENMP
 751                    dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
 752                    dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
 753 #endif
 754                  end if
 755                end do
 756              end do
 757 #if defined HAVE_OPENMP
 758              if(jlmn<nlmn) then
 759                do ilmn=jlmn+1,nlmn
 760                  i0lmn=ilmn*(ilmn-1)/2
 761                  ijlmn=i0lmn+jlmn
 762                  enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
 763                  if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
 764                  do mu=1,ndgxdtfac
 765                    if(cplex_dgxdt(mu)==2)then
 766                      cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor)
 767                    else
 768                      cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
 769                    end if
 770                    dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
 771                    dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
 772                    if (cplex_==2) then
 773                      dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
 774                      dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
 775                    end if
 776                  end do
 777                end do
 778              end if
 779 #endif
 780            end do
 781 !$OMP END DO
 782          end do
 783        end do
 784        ABI_DEALLOCATE(gxfj)
 785 !$OMP END PARALLEL
 786      end if !nspinortot
 787    end if !complex
 788 
 789 !  === Off-diagonal term(s) (up-down, down-up)
 790 
 791 !  --- No parallelization over spinors ---
 792    if (nspinortot==2.and.nspinor==nspinortot) then
 793      ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
 794      ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
 795 !$OMP PARALLEL &
 796 !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), &
 797 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfi,gxfj,ilmn,ijlmn)
 798      ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac))
 799      do ispinor=1,nspinor
 800        jspinor=3-ispinor
 801        do ia=1,nincat
 802          index_enl=atindx1(iatm+ia)
 803 !$OMP DO
 804          do jlmn=1,nlmn
 805            j0lmn=jlmn*(jlmn-1)/2
 806            jjlmn=j0lmn+jlmn
 807            enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor)
 808            do mu=1,ndgxdtfac
 809              if(cplex_dgxdt(mu)==2)then
 810                cplex_ = 2 ;
 811                gxfi(1)    = zero ; gxfi(2)    = dgxdt(1,mu,jlmn,ia,ispinor)
 812                gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,jspinor)
 813              else
 814                cplex_ = cplex ;
 815                gxfi(1:cplex)   =dgxdt(1:cplex,mu,jlmn,ia,ispinor)
 816                gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,jspinor)
 817              end if
 818              dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
 819              dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
 820              if (cplex_==2) then
 821                dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
 822                dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
 823              end if
 824            end do
 825            do ilmn=1,jlmn-1
 826              ijlmn=j0lmn+ilmn
 827              enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
 828              do mu=1,ndgxdtfac
 829                if(cplex_dgxdt(mu)==2)then
 830                  cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor)
 831                else
 832                  cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
 833                end if
 834                dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
 835                dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
 836 #if !defined HAVE_OPENMP
 837                dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
 838                dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
 839 #endif
 840                if (cplex_==2) then
 841                  dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
 842                  dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
 843 #if !defined HAVE_OPENMP
 844                  dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
 845                  dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
 846 #endif
 847                end if
 848              end do !mu
 849            end do !ilmn
 850 #if defined HAVE_OPENMP
 851            if(jlmn<nlmn) then
 852              do ilmn=jlmn+1,nlmn
 853                i0lmn=ilmn*(ilmn-1)/2
 854                ijlmn=i0lmn+jlmn
 855                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
 856                do mu=1,ndgxdtfac
 857                  if(cplex_dgxdt(mu)==2)then
 858                    cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,jspinor)
 859                  else
 860                    cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,jspinor)
 861                  end if
 862                  dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
 863                  dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
 864                  if (cplex_==2) then
 865                    dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
 866                    dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
 867                  end if
 868                end do !mu
 869              end do !ilmn
 870            end if
 871 #endif
 872          end do !jmln
 873 !$OMP END DO
 874        end do !ia
 875      end do !ispinor
 876      ABI_DEALLOCATE(gxfj)
 877 !$OMP END PARALLEL
 878 
 879 !    --- Parallelization over spinors ---
 880    else if (nspinortot==2.and.nspinor/=nspinortot) then
 881      ABI_CHECK(cplex_enl==2,"BUG: opernlc_ylm: invalid cplex_enl/=2!")
 882      ABI_CHECK(cplex_fac==2,"BUG: opernlc_ylm: invalid cplex_fac/=2!")
 883      ABI_ALLOCATE(dgxdtfac_offdiag,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinortot))
 884 !$OMP PARALLEL &
 885 !$OMP PRIVATE(ia,index_enl), &
 886 !$OMP PRIVATE(jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,mu,gxfi)
 887 !$OMP WORKSHARE
 888      dgxdtfac_offdiag(:,:,:,:,:)=zero
 889 !$OMP END WORKSHARE
 890 !$OMP SINGLE
 891      ispinor_index=mpi_enreg%me_spinor+1
 892      jspinor_index=3-ispinor_index
 893      if (ispinor_index==1) then
 894        ijspin=3;jispin=4
 895      else
 896        ijspin=4;jispin=3
 897      end if
 898 !$OMP END SINGLE
 899      do ia=1,nincat
 900        index_enl=atindx1(iatm+ia)
 901 !$OMP DO
 902        do jlmn=1,nlmn
 903          j0lmn=jlmn*(jlmn-1)/2
 904          do ilmn=1,nlmn
 905            i0lmn=ilmn*(ilmn-1)/2
 906            if (ilmn<=jlmn) then
 907              ijlmn=j0lmn+ilmn
 908              enl_(1)= enl_ptr(2*ijlmn-1,index_enl,ijspin)
 909              enl_(2)=-enl_ptr(2*ijlmn  ,index_enl,ijspin)
 910            else
 911              jilmn=i0lmn+jlmn
 912              enl_(1)= enl_ptr(2*jilmn-1,index_enl,jispin)
 913              enl_(2)= enl_ptr(2*jilmn  ,index_enl,jispin)
 914            end if
 915            do mu=1,ndgxdtfac
 916              if(cplex_dgxdt(mu)==2)then
 917                cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
 918              else
 919                cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
 920              end if
 921              dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
 922 &                 dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(1)
 923              dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
 924 &                 dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(2)*gxfi(1)
 925              if (cplex_==2) then
 926                dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
 927 &                   dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)-enl_(2)*gxfi(2)
 928                dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
 929 &                   dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(2)
 930              end if
 931            end do
 932          end do !ilmn
 933        end do !jlmn
 934 !$OMP END DO
 935      end do !iat
 936 !$OMP SINGLE
 937      call xmpi_sum(dgxdtfac_offdiag,mpi_enreg%comm_spinor,ierr)
 938 !$OMP END SINGLE
 939 !$OMP WORKSHARE
 940      dgxdtfac_(:,:,:,:,1)=dgxdtfac_(:,:,:,:,1)+dgxdtfac_offdiag(:,:,:,:,ispinor_index)
 941 !$OMP END WORKSHARE
 942 !$OMP END PARALLEL
 943      ABI_DEALLOCATE(dgxdtfac_offdiag)
 944    end if !nspinortot
 945 
 946   end if ! pawopt & optder
 947 
 948 !Accumulate d2gxdtfac related to nonlocal operator (Norm-conserving)
 949 !-------------------------------------------------------------------
 950   if (optder==2.and.paw_opt==0) then
 951    !Enl is E(Kleinman-Bylander)
 952 !$OMP PARALLEL &
 953 !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_,mu)
 954    do ispinor=1,nspinor
 955      ispinor_index = ispinor + shift
 956      do ia=1,nincat
 957 !$OMP DO
 958        do ilmn=1,nlmn
 959          iln=indlmn(5,ilmn)
 960          enl_(1)=enl_ptr(iln,itypat,ispinor_index)
 961          do mu=1,nd2gxdtfac
 962            d2gxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=enl_(1)*d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
 963          end do
 964        end do
 965 !$OMP END DO
 966      end do
 967    end do
 968 !$OMP END PARALLEL
 969  end if
 970 
 971  DBG_EXIT("COLL")
 972 
 973 !Accumulate d2gxdtfac related to nonlocal operator (PAW)
 974 !-------------------------------------------------------------------
 975   if (optder==2.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then
 976    !Enl is psp strength Dij or (Dij-lambda.Sij)
 977 
 978 !  === Diagonal term(s) (up-up, down-down)
 979 
 980 !  1-Enl is real
 981    if (cplex_enl==1) then
 982 !$OMP PARALLEL &
 983 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi)
 984      ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac))
 985      do ispinor=1,nspinor
 986        ispinor_index=ispinor+shift
 987        do ia=1,nincat
 988          index_enl=atindx1(iatm+ia)
 989 !$OMP DO
 990          do jlmn=1,nlmn
 991            j0lmn=jlmn*(jlmn-1)/2
 992            jjlmn=j0lmn+jlmn
 993            enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
 994            if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
 995            do mu=1,nd2gxdtfac
 996              gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
 997              d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
 998            end do
 999            do ilmn=1,jlmn-1
1000              ijlmn=j0lmn+ilmn
1001              enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
1002              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1003              do mu=1,nd2gxdtfac
1004                gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1005                d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
1006 #if !defined HAVE_OPENMP
1007                d2gxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
1008 #endif
1009              end do
1010            end do
1011 #if defined HAVE_OPENMP
1012            if(jlmn<nlmn) then
1013              do ilmn=jlmn+1,nlmn
1014                i0lmn=ilmn*(ilmn-1)/2
1015                ijlmn=i0lmn+jlmn
1016                enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
1017                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1018                do mu=1,nd2gxdtfac
1019                  gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1020                  d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
1021                end do
1022              end do
1023            end if
1024 #endif
1025          end do
1026 !$OMP END DO
1027        end do
1028      end do
1029      ABI_DEALLOCATE(gxfj)
1030 !$OMP END PARALLEL
1031 
1032 !    2-Enl is complex  ===== D^ss'_ij=D^s's_ji^*
1033    else
1034      ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
1035 
1036      if (nspinortot==1) then ! -------------> NO SPINORS
1037 
1038 !$OMP PARALLEL &
1039 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi)
1040        ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac))
1041        do ia=1,nincat
1042          index_enl=atindx1(iatm+ia)
1043 !$OMP DO
1044          do jlmn=1,nlmn
1045            j0lmn=jlmn*(jlmn-1)/2
1046            jjlmn=j0lmn+jlmn
1047            enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
1048            if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
1049            do mu=1,nd2gxdtfac
1050              if(cplex_d2gxdt(mu)==2)then
1051                cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,1)
1052              else
1053                cplex_ = cplex ; gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,1)
1054              end if
1055              d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu)
1056              if (cplex_==2) d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu)
1057            end do
1058            do ilmn=1,jlmn-1
1059              ijlmn=j0lmn+ilmn
1060              enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
1061              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1062              do mu=1,nd2gxdtfac
1063                if(cplex_d2gxdt(mu)==2)then
1064                  cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1)
1065                else
1066                  cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1)
1067                end if
1068                d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
1069                d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)-enl_(2)*gxfi(1)
1070 #if !defined HAVE_OPENMP
1071                d2gxdtfac_(1,mu,ilmn,ia,1)=d2gxdtfac_(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu)
1072                d2gxdtfac_(2,mu,ilmn,ia,1)=d2gxdtfac_(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu)
1073 #endif
1074                if (cplex_==2) then
1075                  d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(2)*gxfi(2)
1076                  d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
1077 #if !defined HAVE_OPENMP
1078                  d2gxdtfac_(1,mu,ilmn,ia,1)=d2gxdtfac_(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu)
1079                  d2gxdtfac_(2,mu,ilmn,ia,1)=d2gxdtfac_(2,mu,ilmn,ia,1)+enl_(1)*gxfj(2,mu)
1080 #endif
1081                end if
1082              end do
1083            end do
1084 #if defined HAVE_OPENMP
1085            if(jlmn<nlmn) then
1086              do ilmn=jlmn+1,nlmn
1087                i0lmn=ilmn*(ilmn-1)/2
1088                ijlmn=i0lmn+jlmn
1089                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
1090                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1091                do mu=1,nd2gxdtfac
1092                  if(cplex_d2gxdt(mu)==2)then
1093                    cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1)
1094                  else
1095                    cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1)
1096                  end if
1097                  d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
1098                  d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1)
1099                  if (cplex_==2) then
1100                    d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2)
1101                    d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
1102                  end if
1103                end do
1104              end do
1105            end if
1106 #endif
1107          end do
1108 !$OMP END DO
1109        end do
1110        ABI_DEALLOCATE(gxfj)
1111 !$OMP END PARALLEL
1112 
1113      else ! -------------> SPINORIAL CASE
1114 
1115 !  === Diagonal term(s) (up-up, down-down)
1116 
1117 !$OMP PARALLEL &
1118 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
1119 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi)
1120        ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac))
1121        do ispinor=1,nspinor
1122          ispinor_index = ispinor + shift
1123          do ia=1,nincat
1124            index_enl=atindx1(iatm+ia)
1125 !$OMP DO
1126            do jlmn=1,nlmn
1127              j0lmn=jlmn*(jlmn-1)/2
1128              jjlmn=j0lmn+jlmn
1129              enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
1130              if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
1131              do mu=1,nd2gxdtfac
1132                if(cplex_d2gxdt(mu)==2)then
1133                  cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,ispinor)
1134                else
1135                  cplex_ = cplex ; gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1136                end if
1137                d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
1138                if (cplex_==2) d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
1139              end do
1140              do ilmn=1,jlmn-1
1141                ijlmn=j0lmn+ilmn
1142                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
1143                if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1144                do mu=1,nd2gxdtfac
1145                  if(cplex_d2gxdt(mu)==2)then
1146                    cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor)
1147                  else
1148                    cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1149                  end if
1150                  d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
1151                  d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(1)
1152 #if !defined HAVE_OPENMP
1153                  d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
1154                  d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
1155 #endif
1156                  if (cplex_==2) then
1157                    d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(2)
1158                    d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
1159 #if !defined HAVE_OPENMP
1160                    d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
1161                    d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
1162 #endif
1163                  end if
1164                end do
1165              end do
1166 #if defined HAVE_OPENMP
1167              if(jlmn<nlmn) then
1168                do ilmn=jlmn+1,nlmn
1169                  i0lmn=ilmn*(ilmn-1)/2
1170                  ijlmn=i0lmn+jlmn
1171                  enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
1172                  if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1173                  do mu=1,nd2gxdtfac
1174                    if(cplex_d2gxdt(mu)==2)then
1175                      cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor)
1176                    else
1177                      cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1178                    end if
1179                    d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
1180                    d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
1181                    if (cplex_==2) then
1182                      d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
1183                      d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
1184                    end if
1185                  end do
1186                end do
1187              end if
1188 #endif
1189            end do
1190 !$OMP END DO
1191          end do
1192        end do
1193        ABI_DEALLOCATE(gxfj)
1194 !$OMP END PARALLEL
1195      end if !nspinortot
1196    end if !complex
1197 
1198 !  === Off-diagonal term(s) (up-down, down-up)
1199 
1200 !  --- No parallelization over spinors ---
1201    if (nspinortot==2.and.nspinor==nspinortot) then
1202      ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
1203      ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
1204 !$OMP PARALLEL &
1205 !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), &
1206 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfi,gxfj,ilmn,ijlmn)
1207      ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac))
1208      do ispinor=1,nspinor
1209        jspinor=3-ispinor
1210        do ia=1,nincat
1211          index_enl=atindx1(iatm+ia)
1212 !$OMP DO
1213          do jlmn=1,nlmn
1214            j0lmn=jlmn*(jlmn-1)/2
1215            jjlmn=j0lmn+jlmn
1216            enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor)
1217            do mu=1,nd2gxdtfac
1218              if(cplex_d2gxdt(mu)==2)then
1219                cplex_ = 2 ;
1220                gxfi(1)    = zero ; gxfi(2)    = d2gxdt(1,mu,jlmn,ia,ispinor)
1221                gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,jspinor)
1222              else
1223                cplex_ = cplex ;
1224                gxfi(1:cplex)   =d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1225                gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,jspinor)
1226              end if
1227              d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
1228              d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
1229              if (cplex_==2) then
1230                d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
1231                d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
1232              end if
1233            end do
1234            do ilmn=1,jlmn-1
1235              ijlmn=j0lmn+ilmn
1236              enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
1237              do mu=1,nd2gxdtfac
1238                if(cplex_d2gxdt(mu)==2)then
1239                  cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor)
1240                else
1241                  cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1242                end if
1243                d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
1244                d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
1245 #if !defined HAVE_OPENMP
1246                d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
1247                d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
1248 #endif
1249                if (cplex_==2) then
1250                  d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
1251                  d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
1252 #if !defined HAVE_OPENMP
1253                  d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
1254                  d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
1255 #endif
1256                end if
1257              end do !mu
1258            end do !ilmn
1259 #if defined HAVE_OPENMP
1260            if(jlmn<nlmn) then
1261              do ilmn=jlmn+1,nlmn
1262                i0lmn=ilmn*(ilmn-1)/2
1263                ijlmn=i0lmn+jlmn
1264                enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
1265                do mu=1,nd2gxdtfac
1266                  if(cplex_d2gxdt(mu)==2)then
1267                    cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,jspinor)
1268                  else
1269                    cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,jspinor)
1270                  end if
1271                  d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
1272                  d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
1273                  if (cplex_==2) then
1274                    d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
1275                    d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
1276                  end if
1277                end do !mu
1278              end do !ilmn
1279            end if
1280 #endif
1281          end do !jmln
1282 !$OMP END DO
1283        end do !ia
1284      end do !ispinor
1285      ABI_DEALLOCATE(gxfj)
1286 !$OMP END PARALLEL
1287 
1288 !    --- Parallelization over spinors ---
1289    else if (nspinortot==2.and.nspinor/=nspinortot) then
1290      ABI_CHECK(cplex_enl==2,"BUG: opernlc_ylm: invalid cplex_enl/=2!")
1291      ABI_CHECK(cplex_fac==2,"BUG: opernlc_ylm: invalid cplex_fac/=2!")
1292      ABI_ALLOCATE(d2gxdtfac_offdiag,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinortot))
1293 !$OMP PARALLEL &
1294 !$OMP PRIVATE(ia,index_enl), &
1295 !$OMP PRIVATE(jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,mu,gxfi)
1296 !$OMP WORKSHARE
1297      d2gxdtfac_offdiag(:,:,:,:,:)=zero
1298 !$OMP END WORKSHARE
1299 !$OMP SINGLE
1300      ispinor_index=mpi_enreg%me_spinor+1
1301      jspinor_index=3-ispinor_index
1302      if (ispinor_index==1) then
1303        ijspin=3;jispin=4
1304      else
1305        ijspin=4;jispin=3
1306      end if
1307 !$OMP END SINGLE
1308      do ia=1,nincat
1309        index_enl=atindx1(iatm+ia)
1310 !$OMP DO
1311        do jlmn=1,nlmn
1312          j0lmn=jlmn*(jlmn-1)/2
1313          do ilmn=1,nlmn
1314            i0lmn=ilmn*(ilmn-1)/2
1315            if (ilmn<=jlmn) then
1316              ijlmn=j0lmn+ilmn
1317              enl_(1)= enl_ptr(2*ijlmn-1,index_enl,ijspin)
1318              enl_(2)=-enl_ptr(2*ijlmn  ,index_enl,ijspin)
1319            else
1320              jilmn=i0lmn+jlmn
1321              enl_(1)= enl_ptr(2*jilmn-1,index_enl,jispin)
1322              enl_(2)= enl_ptr(2*jilmn  ,index_enl,jispin)
1323            end if
1324            do mu=1,nd2gxdtfac
1325              if(cplex_d2gxdt(mu)==2)then
1326                cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1)
1327              else
1328                cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1)
1329              end if
1330              d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
1331 &                  d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(1)
1332              d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
1333 &                  d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(2)*gxfi(1)
1334              if (cplex_==2) then
1335                d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
1336 &                    d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)-enl_(2)*gxfi(2)
1337                d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
1338 &                    d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(2)
1339              end if
1340            end do
1341          end do !ilmn
1342        end do !jlmn
1343 !$OMP END DO
1344      end do !iat
1345 !$OMP SINGLE
1346      call xmpi_sum(d2gxdtfac_offdiag,mpi_enreg%comm_spinor,ierr)
1347 !$OMP END SINGLE
1348 !$OMP WORKSHARE
1349      d2gxdtfac_(:,:,:,:,1)=d2gxdtfac_(:,:,:,:,1)+d2gxdtfac_offdiag(:,:,:,:,ispinor_index)
1350 !$OMP END WORKSHARE
1351 !$OMP END PARALLEL
1352      ABI_DEALLOCATE(d2gxdtfac_offdiag)
1353    end if !nspinortot
1354 
1355   end if ! pawopt & optder
1356 
1357 !End of loop when a exp(-iqR) phase is present
1358 !------------------------------------------- ------------------------
1359 
1360 !When iphase=1, gxfac and gxfac_ point to the same memory space
1361 !When iphase=2, we add i.gxfac_ to gxfac
1362   if (iphase==2) then
1363 !$OMP PARALLEL PRIVATE(ispinor,ia,ilmn,mu)
1364 !$OMP DO COLLAPSE(3)
1365     do ispinor=1,nspinor
1366       do ia=1,nincat
1367         do ilmn=1,nlmn
1368           gxfac(1,ilmn,ia,ispinor)=gxfac(1,ilmn,ia,ispinor)-gxfac_(2,ilmn,ia,ispinor)
1369           gxfac(2,ilmn,ia,ispinor)=gxfac(2,ilmn,ia,ispinor)+gxfac_(1,ilmn,ia,ispinor)
1370         end do
1371       end do
1372     end do
1373     ABI_DEALLOCATE(gxfac_)
1374     if (optder>=1) then
1375 !$OMP DO COLLAPSE(4)
1376       do ispinor=1,nspinor
1377         do ia=1,nincat
1378           do ilmn=1,nlmn
1379             do mu=1,ndgxdtfac
1380               dgxdtfac(1,mu,ilmn,ia,ispinor)=dgxdtfac(1,mu,ilmn,ia,ispinor)-dgxdtfac_(2,mu,ilmn,ia,ispinor)
1381               dgxdtfac(2,mu,ilmn,ia,ispinor)=dgxdtfac(2,mu,ilmn,ia,ispinor)+dgxdtfac_(1,mu,ilmn,ia,ispinor)
1382             end do
1383           end do
1384         end do
1385       end do
1386       ABI_DEALLOCATE(dgxdtfac_)
1387     end if
1388     if (optder>=2) then
1389 !$OMP DO COLLAPSE(4)
1390       do ispinor=1,nspinor
1391         do ia=1,nincat
1392           do ilmn=1,nlmn
1393             do mu=1,ndgxdtfac
1394               d2gxdtfac(1,mu,ilmn,ia,ispinor)=d2gxdtfac(1,mu,ilmn,ia,ispinor)-d2gxdtfac_(2,mu,ilmn,ia,ispinor)
1395               d2gxdtfac(2,mu,ilmn,ia,ispinor)=d2gxdtfac(2,mu,ilmn,ia,ispinor)+d2gxdtfac_(1,mu,ilmn,ia,ispinor)
1396             end do
1397           end do
1398         end do
1399       end do
1400       ABI_DEALLOCATE(d2gxdtfac_)
1401     end if
1402 !$OMP END PARALLEL
1403   end if
1404 
1405 !End loop over real/imaginary part of the exp(-iqR) phase
1406  end do
1407           
1408 
1409 !Accumulate gxfac related to overlap (Sij) (PAW)
1410 !------------------------------------------- ------------------------
1411  if (paw_opt==3.or.paw_opt==4) then ! Use Sij, overlap contribution
1412 !$OMP PARALLEL &
1413 !$OMP PRIVATE(ispinor,ia,jlmn,j0lmn,jjlmn,jlm,sijr,ilmn,ilm,ijlmn,gxi,gxj)
1414 !$OMP WORKSHARE
1415    gxfac_sij(1:cplex,1:nlmn,1:nincat,1:nspinor)=zero
1416 !$OMP END WORKSHARE
1417 !$OMP DO COLLAPSE(3)
1418    do ispinor=1,nspinor
1419      do ia=1,nincat
1420        do jlmn=1,nlmn
1421          j0lmn=jlmn*(jlmn-1)/2
1422          jjlmn=j0lmn+jlmn
1423          jlm=indlmn(4,jlmn)
1424          sijr=sij(jjlmn)
1425          gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
1426          gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxj(1:cplex)
1427          do ilmn=1,jlmn-1
1428            ilm=indlmn(4,ilmn)
1429           !if (ilm==jlm) then
1430            ijlmn=j0lmn+ilmn
1431            sijr=sij(ijlmn)
1432            gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
1433            gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex)
1434 #if !defined HAVE_OPENMP
1435            gxfac_sij(1:cplex,ilmn,ia,ispinor)=gxfac_sij(1:cplex,ilmn,ia,ispinor)+sijr*gxj(1:cplex)
1436 #endif
1437           !end if
1438          end do
1439 #if defined HAVE_OPENMP
1440          if(jlmn<nlmn) then
1441            do ilmn=jlmn+1,nlmn
1442              ilm=indlmn(4,ilmn)
1443              !if (ilm==jlm) then
1444              i0lmn=ilmn*(ilmn-1)/2
1445              ijlmn=i0lmn+jlmn
1446              sijr=sij(ijlmn)
1447              gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
1448              gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex)
1449              !end if
1450            end do
1451          end if
1452 #endif
1453        end do
1454      end do
1455    end do
1456 !$OMP END DO
1457 !$OMP END PARALLEL
1458  end if
1459 
1460 !Accumulate dgxdtfac related to overlap (Sij) (PAW)
1461 !-------------------------------------------------------------------
1462  if (optder>=1.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution
1463 !$OMP PARALLEL &
1464 !$OMP PRIVATE(ispinor,ia), &
1465 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,sijr,mu,gxfj,ilmn,ijlmn,gxfi)
1466    ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac))
1467 !$OMP WORKSHARE
1468    dgxdtfac_sij(1:cplex,1:ndgxdtfac,1:nlmn,1:nincat,1:nspinor)=zero
1469 !$OMP END WORKSHARE
1470    do ispinor=1,nspinor
1471      do ia=1,nincat
1472 !$OMP DO
1473        do jlmn=1,nlmn
1474          j0lmn=jlmn*(jlmn-1)/2
1475          jjlmn=j0lmn+jlmn
1476          sijr=sij(jjlmn)
1477          do mu=1,ndgxdtfac
1478            gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
1479            dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1480          end do
1481          do ilmn=1,jlmn-1
1482            ijlmn=j0lmn+ilmn
1483            sijr=sij(ijlmn)
1484            do mu=1,ndgxdtfac
1485              gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
1486              dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1487 #if !defined HAVE_OPENMP
1488              dgxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1489 #endif
1490            end do
1491          end do
1492 #if defined HAVE_OPENMP
1493          if(jlmn<nlmn) then
1494            do ilmn=jlmn+1,nlmn
1495              i0lmn=ilmn*(ilmn-1)/2
1496              ijlmn=i0lmn+jlmn
1497              sijr=sij(ijlmn)
1498              do mu=1,ndgxdtfac
1499                gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
1500                dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1501              end do
1502            end do
1503          end if
1504 #endif
1505        end do
1506 !$OMP END DO
1507      end do
1508    end do
1509    ABI_DEALLOCATE(gxfj)
1510 !$OMP END PARALLEL
1511  end if
1512 
1513 !Accumulate d2gxdtfac related to overlap (Sij) (PAW)
1514 !-------------------------------------------------------------------
1515  if (optder==2.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution
1516 !$OMP PARALLEL &
1517 !$OMP PRIVATE(ispinor,ia), &
1518 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,sijr,mu,gxfj,ilmn,ijlmn,gxfi)
1519    ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac))
1520 !$OMP WORKSHARE
1521    d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,1:nlmn,1:nincat,1:nspinor)=zero
1522 !$OMP END WORKSHARE
1523    do ispinor=1,nspinor
1524      do ia=1,nincat
1525 !$OMP DO
1526        do jlmn=1,nlmn
1527          j0lmn=jlmn*(jlmn-1)/2
1528          jjlmn=j0lmn+jlmn
1529          sijr=sij(jjlmn)
1530          do mu=1,nd2gxdtfac
1531            gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1532            d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)= &
1533 &                d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1534          end do
1535          do ilmn=1,jlmn-1
1536            ijlmn=j0lmn+ilmn
1537            sijr=sij(ijlmn)
1538            do mu=1,nd2gxdtfac
1539              gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1540              d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)= &
1541 &                  d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1542 #if !defined HAVE_OPENMP
1543              d2gxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)= &
1544 &                  d2gxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1545 #endif
1546            end do
1547          end do
1548 #if defined HAVE_OPENMP
1549          if(jlmn<nlmn) then
1550            do ilmn=jlmn+1,nlmn
1551              i0lmn=ilmn*(ilmn-1)/2
1552              ijlmn=i0lmn+jlmn
1553              sijr=sij(ijlmn)
1554              do mu=1,nd2gxdtfac
1555                gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1556                d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)= &
1557 &                    d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1558              end do
1559            end do
1560          end if
1561 #endif
1562        end do
1563 !$OMP END DO
1564      end do
1565    end do
1566    ABI_DEALLOCATE(gxfj)
1567 !$OMP END PARALLEL
1568  end if
1569 
1570 end subroutine opernlc_ylm