TABLE OF CONTENTS


ABINIT/m_opernlb_ylm [ Modules ]

[ Top ] [ Modules ]

NAME

  m_opernlb_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_opernlb_ylm
27 
28  use defs_basis
29  use m_abicore
30  use m_errors
31 #if defined HAVE_OPENMP
32  use OMP_LIB
33 #endif
34 
35  implicit none
36 
37  private

ABINIT/opernlb_ylm [ Functions ]

[ Top ] [ Functions ]

NAME

 opernlb_ylm

FUNCTION

 * Operate with the non-local part of the hamiltonian,
   from projected scalars to reciprocal space.
 * Operate with the non-local projectors and the overlap matrix,
   from projected scalars to reciprocal space.

INPUTS

  choice=chooses possible output (see below)
  cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
        2 if <p_lmn|c> scalars are complex
  cplex_dgxdt(ndgxdt_fac) = used only when cplex = 1
    cplex_dgxdt(i)=1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
  cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
  dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Vnl (NL operator)
  dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfacrelated to Sij (overlap)
  dimffnl=second dimension of ffnl
  ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
  gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
  gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))= reduced projected scalars related to Sij (overlap)
  ia3=gives the number of the first atom in the subset presently treated
  idir=direction of the - atom to be moved in the case (choice=2,signs=2),
                        - k point direction in the case (choice=5, 51, 52 and signs=2)
                        - strain component (1:6) in the case (choice=2,signs=2) or (choice=6,signs=1)
  indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
  kpg(npw,nkpg)=(k+G) components (if nkpg=3)
  matblk=dimension of the array ph3d
  ndgxdtfac=second dimension of dgxdtfac
  nincat=number of atoms in the subset here treated
  nkpg=second dimension of array kpg (0 or 3)
  nlmn=number of (l,m,n) numbers for current type of atom
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npw=number of plane waves in reciprocal space
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  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)
  ph3d(2,npw,matblk)=three-dimensional phase factors
  ucvol=unit cell volume (bohr^3)

OUTPUT

  (see side effects)

SIDE EFFECTS

 --if (paw_opt=0, 1 or 4)
    vect(2,npwout*nspinor)=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=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|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)
    vect(2,npwout*nspinor)=final vector in reciprocal space:
      if (choice=1)  <G|V_nonlocal-lamdba.(I+S)|vect_in> (note: not including <G|I|c>)
      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=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|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)
    svect(2,npwout*nspinor)=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> (note: not including <G|I|c>)
      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=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|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>

NOTES

 1-The openMP version is different from the standard version:
   the standard version is more effifient on one CPU core.
 2-Operate for one type of atom, and within this given type of atom,
   for a subset of at most nincat atoms.

PARENTS

      nonlop_ylm

CHILDREN

SOURCE

 150 subroutine opernlb_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,&
 151 &                      d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnl,ffnl,gxfac,gxfac_sij,&
 152 &                      ia3,idir,indlmn,kpg,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpg,nlmn,nloalg,npw,&
 153 &                      nspinor,paw_opt,ph3d,svect,ucvol,vect)
 154 
 155 
 156 !This section has been created automatically by the script Abilint (TD).
 157 !Do not modify the following lines by hand.
 158 #undef ABI_FUNC
 159 #define ABI_FUNC 'opernlb_ylm'
 160 !End of the abilint section
 161 
 162  implicit none
 163 
 164 !Arguments ------------------------------------
 165 !scalars
 166 
 167  integer,intent(in) :: choice,cplex,cplex_fac,dimffnl,ia3,idir,matblk,ndgxdtfac,nd2gxdtfac,nincat
 168  integer,intent(in) :: nkpg,nlmn,npw,nspinor,paw_opt
 169  real(dp),intent(in) :: ucvol
 170 !arrays
 171  integer,intent(in) ::  cplex_dgxdt(ndgxdtfac),cplex_d2gxdt(nd2gxdtfac),indlmn(6,nlmn),nloalg(3)
 172  real(dp),intent(in) :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor)
 173  real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)
 174  real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat*(paw_opt/3),nspinor)
 175  real(dp),intent(in) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat*(paw_opt/3),nspinor)
 176  real(dp),intent(in) :: ffnl(npw,dimffnl,nlmn),gxfac(cplex_fac,nlmn,nincat,nspinor)
 177  real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
 178  real(dp),intent(in) :: kpg(npw,nkpg),ph3d(2,npw,matblk)
 179  real(dp),intent(inout) :: svect(:,:),vect(:,:)
 180 !Local variables-------------------------------
 181 !Arrays
 182 !scalars
 183  integer :: fdb,fdf,ia,iaph3d,ic,ii,il,ilmn,ipw,ipwshft,ispinor,jc,nthreads,ffnl_dir1,ffnl_dir(3)
 184  real(dp) :: scale,wt
 185  logical :: parity
 186 !arrays
 187  integer,parameter :: ffnl_dir_dat(6)=(/3,4,4,2,2,3/)
 188  integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
 189  integer,parameter :: idir1(9)=(/1,1,1,2,2,2,3,3,3/),idir2(9)=(/1,2,3,1,2,3,1,2,3/)
 190  real(dp),allocatable :: d2gxdtfac_(:,:,:),d2gxdtfacs_(:,:,:),dgxdtfac_(:,:,:),dgxdtfacs_(:,:,:),gxfac_(:,:),gxfacs_(:,:)
 191  complex(dpc),allocatable :: ztab(:)
 192 
 193 ! *************************************************************************
 194 
 195  DBG_ENTER("COLL")
 196 
 197 !Nothing to do when choice=4, 6 or 23
 198  if (choice==4.or.choice==6.or.choice==23) return
 199 
 200 !DDK not compatible with istwkf > 1
 201  if(cplex==1.and.(any(cplex_dgxdt(:)==2).or.any(cplex_d2gxdt(:)==2)))then
 202    MSG_BUG("opernlb_ylm+ddk not compatible with istwfk>1")
 203  end if
 204 
 205 !Inits
 206  wt=four_pi/sqrt(ucvol)
 207  nthreads=1
 208 #if defined HAVE_OPENMP
 209  nthreads=OMP_GET_NUM_THREADS()
 210 #endif
 211 
 212  if (paw_opt/=3) then
 213    ABI_ALLOCATE(gxfac_,(2,nlmn))
 214    gxfac_(:,:)=zero
 215    if (choice>1) then
 216      ABI_ALLOCATE(dgxdtfac_,(2,ndgxdtfac,nlmn))
 217      if(ndgxdtfac>0) dgxdtfac_(:,:,:)=zero
 218    end if
 219    if (choice==54.or.choice==8.or.choice==81) then
 220      ABI_ALLOCATE(d2gxdtfac_,(2,nd2gxdtfac,nlmn))
 221      if(nd2gxdtfac>0) d2gxdtfac_(:,:,:)=zero
 222    end if
 223  end if
 224  if (paw_opt>=3) then
 225    ABI_ALLOCATE(gxfacs_,(2,nlmn))
 226    gxfacs_(:,:)=zero
 227    if (choice>1) then
 228      ABI_ALLOCATE(dgxdtfacs_,(2,ndgxdtfac,nlmn))
 229      if (ndgxdtfac>0) dgxdtfacs_(:,:,:)=zero
 230    end if
 231    if (choice==54.or.choice==8.or.choice==81) then
 232      ABI_ALLOCATE(d2gxdtfacs_,(2,nd2gxdtfac,nlmn))
 233      if (nd2gxdtfac>0) d2gxdtfacs_(:,:,:)=zero
 234    end if
 235  end if
 236 
 237  ABI_ALLOCATE(ztab,(npw))
 238 
 239 !==========================================================================
 240 !========== STANDARD VERSION ==============================================
 241 !==========================================================================
 242  if (nthreads==1) then
 243 
 244 !  Loop on spinorial components
 245    do ispinor=1,nspinor
 246      ipwshft=(ispinor-1)*npw
 247 
 248 !    Loop on atoms (blocking)
 249      do ia=1,nincat
 250        iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
 251 !      Scale gxfac with 4pi/sqr(omega).(-i)^l
 252        if (paw_opt/=3) then
 253          do ilmn=1,nlmn
 254            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 255            scale=wt;if (il>1) scale=-scale
 256            if (parity) then
 257              gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor)
 258              if (cplex_fac==1) gxfac_(2,ilmn)=zero
 259            else
 260              gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor)
 261              if (cplex_fac==2) then
 262                gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor)
 263              else
 264                gxfac_(1,ilmn)=zero
 265              end if
 266            end if
 267          end do
 268          if (choice>1) then
 269            do ilmn=1,nlmn
 270              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 271              scale=wt;if (il>1) scale=-scale
 272              if (parity) then
 273                if(cplex_fac==2)then
 274                  dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor)
 275                else
 276                  do ii=1,ndgxdtfac
 277                    ic = cplex_dgxdt(ii) ; jc = 3-ic
 278                    dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 279                    dgxdtfac_(jc,ii,ilmn)=zero
 280                  end do
 281                end if
 282              else
 283                if(cplex_fac==2)then
 284                  do ii=1,ndgxdtfac
 285                    dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor)
 286                    dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 287                  end do
 288                else
 289                  do ii=1,ndgxdtfac
 290                    ic =  cplex_dgxdt(ii) ; jc = 3-ic
 291                    dgxdtfac_(ic,ii,ilmn)=zero
 292                    if(ic==1)then
 293                      dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 294                    else
 295                      dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 296                    end if
 297                  end do
 298                end if
 299              end if
 300            end do
 301          end if
 302          if (choice==54.or.choice==8.or.choice==81) then
 303            do ilmn=1,nlmn
 304              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 305              scale=wt;if (il>1) scale=-scale
 306              if (parity) then
 307                if(cplex_fac==2)then
 308                  d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor)
 309                else
 310                  do ii=1,nd2gxdtfac
 311                    ic = cplex_d2gxdt(ii) ; jc = 3-ic
 312                    d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 313                    d2gxdtfac_(jc,ii,ilmn)=zero
 314                  end do
 315                end if
 316              else
 317                if(cplex_fac==2)then
 318                  do ii=1,nd2gxdtfac
 319                    d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor)
 320                    d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 321                  end do
 322                else
 323                  do ii=1,nd2gxdtfac
 324                    ic =  cplex_d2gxdt(ii) ; jc = 3-ic
 325                    d2gxdtfac_(ic,ii,ilmn)=zero
 326                    if(ic==1)then
 327                      d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 328                    else
 329                      d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 330                    end if
 331                  end do
 332                end if
 333              end if
 334            end do
 335          end if
 336        end if
 337 
 338 !      Scale gxfac_sij with 4pi/sqr(omega).(-i)^l
 339        if (paw_opt>=3) then
 340          do ilmn=1,nlmn
 341            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 342            scale=wt;if (il>1) scale=-scale
 343            if (parity) then
 344              gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor)
 345              if (cplex==1) gxfacs_(2,ilmn)=zero
 346            else
 347              gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor)
 348              if (cplex==2) then
 349                gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor)
 350              else
 351                gxfacs_(1,ilmn)=zero
 352              end if
 353            end if
 354          end do
 355          if (choice>1) then
 356            do ilmn=1,nlmn
 357              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 358              scale=wt;if (il>1) scale=-scale
 359              if (parity) then
 360                if(cplex==2)then
 361                  dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn)=scale*dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor)
 362                else
 363                  do ii=1,ndgxdtfac
 364                    ic = cplex_dgxdt(ii) ; jc = 3-ic
 365                    dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 366                    dgxdtfacs_(jc,ii,ilmn)=zero
 367                  end do
 368                end if
 369              else
 370                if(cplex==2)then
 371                  do ii=1,ndgxdtfac
 372                    dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor)
 373                    dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 374                  end do
 375                else
 376                  do ii=1,ndgxdtfac
 377                    ic = cplex_dgxdt(ii) ; jc = 3-ic
 378                    dgxdtfacs_(ic,ii,ilmn)=zero
 379                    if(ic==1)then
 380                      dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 381                    else
 382                      dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 383                    end if
 384                  end do
 385                end if
 386              end if
 387            end do
 388          end if
 389          if (choice==54.or.choice==8.or.choice==81) then
 390            do ilmn=1,nlmn
 391              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 392              scale=wt;if (il>1) scale=-scale
 393              if (parity) then
 394                if(cplex==2)then
 395                  d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor)
 396                else
 397                  do ii=1,nd2gxdtfac
 398                    ic = cplex_d2gxdt(ii) ; jc = 3-ic
 399                    d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 400                    d2gxdtfacs_(jc,ii,ilmn)=zero
 401                  end do
 402                end if
 403              else
 404                if(cplex==2)then
 405                  do ii=1,nd2gxdtfac
 406                    d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor)
 407                    d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 408                  end do
 409                else
 410                  do ii=1,nd2gxdtfac
 411                    ic = cplex_d2gxdt(ii) ; jc = 3-ic
 412                    d2gxdtfacs_(ic,ii,ilmn)=zero
 413                    if(ic==1)then
 414                      d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 415                    else
 416                      d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 417                    end if
 418                  end do
 419                end if
 420              end if
 421            end do
 422          end if
 423        end if
 424 
 425 !      Compute <g|Vnl|c> (or derivatives) for each plane wave:
 426 
 427        if (paw_opt/=3) then
 428 
 429          ztab(:)=czero
 430 
 431 !        ------
 432          if (choice==1) then ! <g|Vnl|c>
 433            do ilmn=1,nlmn
 434              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 435            end do
 436          end if
 437 
 438 !        ------
 439          if (choice==2) then ! derivative w.r.t. atm. pos
 440            do ilmn=1,nlmn
 441              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp)
 442            end do
 443            ztab(:)=two_pi*kpg(:,idir)*ztab(:)
 444            do ilmn=1,nlmn
 445              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
 446            end do
 447          end if
 448 
 449 !        ------
 450          if (choice==3) then ! derivative w.r.t. strain
 451            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 452            if (idir<=3) then
 453              do ilmn=1,nlmn
 454                ztab(:)=ztab(:)+ffnl(:,1,ilmn)&
 455 &               *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp)&
 456 &               -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 457              end do
 458            else
 459              do ilmn=1,nlmn
 460                ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)&
 461 &               -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 462              end do
 463            end if
 464          end if
 465 
 466 !        ------
 467          if (choice==5) then ! full derivative w.r.t. k
 468            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 469            do ilmn=1,nlmn
 470              ztab(:)=ztab(:)+ffnl(:,1        ,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)&
 471 &             +ffnl(:,ffnl_dir1,ilmn)*cmplx(   gxfac_(1,  ilmn),   gxfac_(2  ,ilmn),kind=dp)
 472            end do
 473          end if
 474 
 475 !        ------
 476          if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi>
 477            do ilmn=1,nlmn
 478              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
 479            end do
 480          end if
 481 
 482 !        ------
 483          if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi>
 484            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 485            do ilmn=1,nlmn
 486              ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 487            end do
 488          end if
 489 
 490 !        ------
 491          if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>V<dp/dk_(idir-1)|psi>
 492 !                                                -<G|dp/dk_(idir-1)>V<dp/dk_(idir+1)|psi>
 493            fdf = ffnl_dir_dat(2*idir-1)
 494            fdb = ffnl_dir_dat(2*idir)
 495            do ilmn=1,nlmn
 496              ztab(:)=ztab(:) + &
 497 &             ffnl(:,fdf,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) - &
 498 &             ffnl(:,fdb,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
 499            end do
 500          end if
 501 
 502 !        ------
 503          if (choice==54) then ! mixed derivative w.r.t. atm. pos and (right) k
 504            do ilmn=1,nlmn
 505              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp)
 506            end do
 507            ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:)
 508            do ilmn=1,nlmn
 509              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)
 510            end do
 511          end if
 512 
 513 !        ------
 514          if (choice==8) then ! full second order derivative w.r.t. k
 515            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
 516            ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
 517            ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
 518            do ilmn=1,nlmn
 519              ztab(:)=ztab(:) &
 520 &             +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx(    gxfac_(1,  ilmn),    gxfac_(2,  ilmn),kind=dp)&
 521 &             +ffnl(:,1            ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
 522 &             +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)&
 523 &             +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
 524            end do
 525          end if
 526 
 527 !        ------
 528          if (choice==81) then
 529            ! partial second order derivative w.r.t. k
 530            ! full derivative w.r.t. k1, right derivative w.r.t. k2
 531            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
 532            ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
 533            do ilmn=1,nlmn
 534              ztab(:)=ztab(:) &
 535 &             +ffnl(:,1            ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
 536 &             +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
 537            end do
 538          end if
 539 
 540 !        ------
 541          ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
 542 
 543          vect(1,1+ipwshft:npw+ipwshft)=vect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
 544          vect(2,1+ipwshft:npw+ipwshft)=vect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:))
 545 
 546        end if
 547 
 548 !      Compute <g|S|c> (or derivatives) for each plane wave:
 549 
 550        if (paw_opt>=3) then
 551 
 552          ztab(:)=czero
 553 
 554 !        ------
 555          if (choice==1) then ! <g|S|c>
 556            do ilmn=1,nlmn
 557              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
 558            end do
 559          end if
 560 
 561 !        ------
 562          if (choice==2) then ! derivative w.r.t. atm. pos
 563            do ilmn=1,nlmn
 564              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp)
 565            end do
 566            ztab(:)=two_pi*kpg(:,idir)*ztab(:)
 567            do ilmn=1,nlmn
 568              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
 569            end do
 570          end if
 571 
 572 !        ------
 573          if (choice==3) then ! derivative w.r.t. strain
 574            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 575            if (idir<=3) then
 576              do ilmn=1,nlmn
 577                ztab(:)=ztab(:)+ffnl(:,1,ilmn)&
 578 &               *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)&
 579 &               -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
 580              end do
 581            else
 582              do ilmn=1,nlmn
 583                ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)&
 584 &               -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
 585              end do
 586            end if
 587          end if
 588 
 589 !        ------
 590          if (choice==5) then ! full derivative w.r.t. k
 591            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 592            do ilmn=1,nlmn
 593              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)&
 594 &             +ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
 595            end do
 596          end if
 597 
 598 !        ------
 599          if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi>
 600            do ilmn=1,nlmn
 601              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
 602            end do
 603          end if
 604 
 605 !        ------
 606          if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi>
 607            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 608            do ilmn=1,nlmn
 609              ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
 610            end do
 611          end if
 612 
 613 !        ------
 614          if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>S<dp/dk_(idir-1)|psi>
 615 !                                                -<G|dp/dk_(idir-1)>S<dp/dk_(idir+1)|psi>
 616            fdf = ffnl_dir_dat(2*idir-1)
 617            fdb = ffnl_dir_dat(2*idir)
 618            do ilmn=1,nlmn
 619              ztab(:)=ztab(:) + &
 620 &             ffnl(:,fdf,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp) - &
 621 &             ffnl(:,fdb,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
 622            end do
 623          end if
 624 
 625 !        ------
 626          if (choice==54) then ! mixed derivative w.r.t. atm. pos and k
 627            do ilmn=1,nlmn
 628              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp)
 629            end do
 630            ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:)
 631            do ilmn=1,nlmn
 632              ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)
 633            end do
 634          end if
 635 
 636 !        ------
 637          if (choice==8) then ! full second order derivative w.r.t. k
 638            ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
 639            ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
 640            do ilmn=1,nlmn
 641              ztab(:)=ztab(:) &
 642 &             +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx(    gxfacs_(1,  ilmn),    gxfacs_(2,  ilmn),kind=dp)&
 643 &             +ffnl(:,1            ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
 644 &             +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)&
 645 &             +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
 646            end do
 647          end if
 648 
 649 !        ------
 650          if (choice==81) then
 651            ! partial second order derivative w.r.t. k
 652            ! full derivative w.r.t. k1, right derivative w.r.t. k2
 653            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
 654            ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
 655            do ilmn=1,nlmn
 656              ztab(:)=ztab(:) &
 657 &             +ffnl(:,1            ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
 658 &             +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
 659            end do
 660          end if
 661 
 662 
 663 !        ------
 664          ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
 665          svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
 666          svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:))
 667        end if
 668 
 669 !      End loop on atoms
 670      end do
 671    end do !  End loop on spinors
 672 
 673 
 674 !  ==========================================================================
 675 !  ========== OPENMP VERSION ================================================
 676 !  ==========================================================================
 677  else
 678 
 679 !  Loop on spinorial components
 680    do ispinor=1,nspinor
 681      ipwshft=(ispinor-1)*npw
 682 
 683 !    Loop on atoms (blocking)
 684      do ia=1,nincat
 685        iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
 686 
 687 !      Scale gxfac with 4pi/sqr(omega).(-i)^l
 688        if (paw_opt/=3) then
 689 !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc)
 690 !$OMP DO
 691          do ilmn=1,nlmn
 692            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 693            scale=wt;if (il>1) scale=-scale
 694            if (parity) then
 695              gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor)
 696              if (cplex_fac==1) gxfac_(2,ilmn)=zero
 697            else
 698              gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor)
 699              if (cplex_fac==2) then
 700                gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor)
 701              else
 702                gxfac_(1,ilmn)=zero
 703              end if
 704            end if
 705          end do
 706 !$OMP END DO
 707          if (choice>1) then
 708 !$OMP DO
 709            do ilmn=1,nlmn
 710              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 711              scale=wt;if (il>1) scale=-scale
 712              if (parity) then
 713                if(cplex_fac==2)then
 714                  dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor)
 715                else
 716                  do ii=1,ndgxdtfac
 717                    ic = cplex_dgxdt(ii) ; jc = 3-ic
 718                    dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 719                    dgxdtfac_(jc,ii,ilmn)=zero
 720                  end do
 721                end if
 722              else
 723                if(cplex_fac==2)then
 724                  do ii=1,ndgxdtfac
 725                    dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 726                    dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor)
 727                  end do
 728                else
 729                  do ii=1,ndgxdtfac
 730                    ic =  cplex_dgxdt(ii) ; jc = 3-ic
 731                    dgxdtfac_(ic,ii,ilmn)=zero
 732                    if(ic==1)then
 733                      dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 734                    else
 735                      dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
 736                    end if
 737                  end do
 738                end if
 739              end if
 740            end do
 741 !$OMP END DO
 742          end if
 743          if (choice==54.or.choice==8.or.choice==81) then
 744 !$OMP DO
 745            do ilmn=1,nlmn
 746              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 747              scale=wt;if (il>1) scale=-scale
 748              if (parity) then
 749                if(cplex_fac==2)then
 750                  d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor)
 751                else
 752                  do ii=1,nd2gxdtfac
 753                    ic = cplex_d2gxdt(ii) ; jc = 3-ic
 754                    d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 755                    d2gxdtfac_(jc,ii,ilmn)=zero
 756                  end do
 757                end if
 758              else
 759                if(cplex_fac==2)then
 760                  do ii=1,nd2gxdtfac
 761                    d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor)
 762                    d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 763                  end do
 764                else
 765                  do ii=1,nd2gxdtfac
 766                    ic =  cplex_d2gxdt(ii) ; jc = 3-ic
 767                    d2gxdtfac_(ic,ii,ilmn)=zero
 768                    if(ic==1)then
 769                      d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 770                    else
 771                      d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
 772                    end if
 773                  end do
 774                end if
 775              end if
 776            end do
 777 !$OMP END DO
 778          end if
 779 !$OMP END PARALLEL
 780        end if
 781 
 782 !      Scale gxfac_sij with 4pi/sqr(omega).(-i)^l
 783        if (paw_opt>=3) then
 784 !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc)
 785 !$OMP DO
 786          do ilmn=1,nlmn
 787            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 788            scale=wt;if (il>1) scale=-scale
 789            if (parity) then
 790              gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor)
 791              if (cplex==1) gxfacs_(2,ilmn)=zero
 792            else
 793              gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor)
 794              if (cplex==2) then
 795                gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor)
 796              else
 797                gxfacs_(1,ilmn)=zero
 798              end if
 799            end if
 800          end do
 801 !$OMP END DO
 802          if (choice>1) then
 803 !$OMP DO
 804            do ilmn=1,nlmn
 805              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 806              scale=wt;if (il>1) scale=-scale
 807              if (parity) then
 808                if(cplex==2)then
 809                  dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn)=scale*dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor)
 810                else
 811                  do ii=1,ndgxdtfac
 812                    ic = cplex_dgxdt(ii) ; jc = 3-ic
 813                    dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 814                    dgxdtfacs_(jc,ii,ilmn)=zero
 815                  end do
 816                end if
 817              else
 818                if(cplex==2)then
 819                  do ii=1,ndgxdtfac
 820                    dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 821                    dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor)
 822                  end do
 823                else
 824                  do ii=1,ndgxdtfac
 825                    ic = cplex_dgxdt(ii) ; jc = 3-ic
 826                    dgxdtfacs_(ic,ii,ilmn)=zero
 827                    if(ic==1)then
 828                      dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 829                    else
 830                      dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
 831                    end if
 832                  end do
 833                end if
 834              end if
 835            end do
 836 !$OMP END DO
 837          end if
 838          if (choice==54.or.choice==8.or.choice==81) then
 839 !$OMP DO
 840            do ilmn=1,nlmn
 841              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 842              scale=wt;if (il>1) scale=-scale
 843              if (parity) then
 844                if(cplex==2)then
 845                  d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor)
 846                else
 847                  do ii=1,nd2gxdtfac
 848                    ic = cplex_d2gxdt(ii) ; jc = 3-ic
 849                    d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 850                    d2gxdtfacs_(jc,ii,ilmn)=zero
 851                  end do
 852                end if
 853              else
 854                if(cplex==2)then
 855                  do ii=1,nd2gxdtfac
 856                    d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor)
 857                    d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 858                  end do
 859                else
 860                  do ii=1,nd2gxdtfac
 861                    ic = cplex_d2gxdt(ii) ; jc = 3-ic
 862                    d2gxdtfacs_(ic,ii,ilmn)=zero
 863                    if(ic==1)then
 864                      d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 865                    else
 866                      d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
 867                    end if
 868                  end do
 869                end if
 870              end if
 871            end do
 872 !$OMP END DO
 873          end if
 874 !$OMP END PARALLEL
 875        end if
 876 
 877 !      Compute <g|Vnl|c> (or derivatives) for each plane wave:
 878        if (paw_opt/=3) then
 879 !$OMP PARALLEL PRIVATE(ipw,ilmn,fdf,fdb,ffnl_dir1)
 880 
 881 !        ------
 882          if (choice==1) then ! <g|Vnl|c>
 883 !$OMP DO
 884            do ipw=1,npw
 885              ztab(ipw)=czero
 886              do ilmn=1,nlmn
 887                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 888              end do
 889            end do
 890 !$OMP END DO
 891 
 892 !        ------
 893          else if (choice==2) then ! derivative w.r.t. atm. pos
 894 !$OMP DO
 895            do ipw=1,npw
 896              ztab(ipw)=czero
 897              do ilmn=1,nlmn
 898                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp)
 899              end do
 900              ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw)
 901              do ilmn=1,nlmn
 902                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
 903              end do
 904            end do
 905 !$OMP END DO
 906 
 907 !        ------
 908          else if (choice==3) then ! derivative w.r.t. strain
 909            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 910            if (idir<=3) then
 911 !$OMP DO
 912              do ipw=1,npw
 913                ztab(ipw)=czero
 914                do ilmn=1,nlmn
 915                  ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) &
 916 &                 *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp) &
 917 &                 -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 918                end do
 919              end do
 920 !$OMP END DO
 921            else
 922 !$OMP DO
 923              do ipw=1,npw
 924                ztab(ipw)=czero
 925                do ilmn=1,nlmn
 926                  ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) &
 927 &                 -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 928                end do
 929              end do
 930 !$OMP END DO
 931            end if
 932 
 933 
 934 !        ------
 935          else if (choice==5) then ! full derivative w.r.t. k
 936            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 937 !$OMP DO
 938            do ipw=1,npw
 939              ztab(ipw)=czero
 940              do ilmn=1,nlmn
 941                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) &
 942 &               +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 943              end do
 944            end do
 945 !$OMP END DO
 946 
 947 !        ------
 948          else if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi>
 949 !$OMP DO
 950            do ipw=1,npw
 951              ztab(ipw)=czero
 952              do ilmn=1,nlmn
 953                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
 954              end do
 955            end do
 956 !$OMP END DO
 957 
 958 !        ------
 959          else if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi>
 960            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
 961 !$OMP DO
 962            do ipw=1,npw
 963              ztab(ipw)=czero
 964              do ilmn=1,nlmn
 965                ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
 966              end do
 967            end do
 968 !$OMP END DO
 969 
 970 !        ------
 971          else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>V<dp/dk_(idir-1)|psi>
 972 !                                                     -<G|dp/dk_(idir-1)>V<dp/dk_(idir+1)|psi>
 973            fdf = ffnl_dir_dat(2*idir-1)
 974            fdb = ffnl_dir_dat(2*idir)
 975 !$OMP DO
 976            do ipw=1,npw
 977              ztab(ipw)=czero
 978              do ilmn=1,nlmn
 979                ztab(ipw)=ztab(ipw) &
 980 &               +ffnl(ipw,fdf,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) &
 981 &               -ffnl(ipw,fdb,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
 982              end do
 983            end do
 984 !$OMP END DO
 985 
 986 !        ------
 987          else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k
 988 !$OMP DO
 989            do ipw=1,npw
 990              ztab(ipw)=czero
 991              do ilmn=1,nlmn
 992                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp)
 993              end do
 994              ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw)
 995              do ilmn=1,nlmn
 996                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)
 997              end do
 998            end do
 999 !$OMP END DO
1000 
1001 !        ------
1002          else if (choice==8) then ! full second order derivative w.r.t. k
1003            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1004            ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
1005            ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
1006 !$OMP DO
1007            do ipw=1,npw
1008              ztab(ipw)=czero
1009              do ilmn=1,nlmn
1010                ztab(ipw)=ztab(ipw) &
1011 &               +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx(    gxfac_(1,  ilmn),    gxfac_(2,  ilmn),kind=dp)&
1012 &               +ffnl(ipw,1            ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
1013 &               +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)&
1014 &               +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
1015              end do
1016            end do
1017 !$OMP END DO
1018 
1019 !        ------
1020          else if (choice==81) then
1021            ! partial second order derivative w.r.t. k
1022            ! full derivative w.r.t. k1, right derivative w.r.t. k2
1023            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1024            ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
1025 !$OMP DO
1026            do ipw=1,npw
1027              ztab(ipw)=czero
1028              do ilmn=1,nlmn
1029                ztab(ipw)=ztab(ipw) &
1030 &               +ffnl(ipw,1            ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
1031 &               +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
1032              end do
1033            end do
1034 !$OMP END DO
1035 
1036 !        ------
1037          else
1038 !$OMP WORKSHARE
1039            ztab(:)=czero
1040 !$OMP END WORKSHARE
1041          end if
1042 
1043 
1044 !        ------
1045 !$OMP DO
1046          do ipw=1,npw
1047            ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp)
1048            vect(1,ipw+ipwshft)=vect(1,ipw+ipwshft)+real(ztab(ipw))
1049            vect(2,ipw+ipwshft)=vect(2,ipw+ipwshft)+aimag(ztab(ipw))
1050          end do
1051 !$OMP END DO
1052 
1053 !$OMP END PARALLEL
1054        end if
1055 
1056 !      Compute <g|S|c> (or derivatives) for each plane wave:
1057        if (paw_opt>=3) then
1058 !$OMP PARALLEL PRIVATE(ilmn,ipw,fdf,fdb,ffnl_dir1)
1059 
1060 !        ------
1061          if (choice==1) then ! <g|S|c>
1062 !$OMP DO
1063            do ipw=1,npw
1064              ztab(ipw)=czero
1065              do ilmn=1,nlmn
1066                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1067              end do
1068            end do
1069 !$OMP END DO
1070 
1071 !        ------
1072          else if (choice==2) then ! derivative w.r.t. atm. pos
1073 !$OMP DO
1074            do ipw=1,npw
1075              ztab(ipw)=czero
1076              do ilmn=1,nlmn
1077                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp)
1078              end do
1079              ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw)
1080              do ilmn=1,nlmn
1081                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
1082              end do
1083            end do
1084 !$OMP END DO
1085 
1086 !        ------
1087          else if (choice==3) then ! derivative w.r.t. strain
1088            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1089            if (idir<=3) then
1090 !$OMP DO
1091              do ipw=1,npw
1092                ztab(ipw)=czero
1093                do ilmn=1,nlmn
1094                  ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) &
1095 &                 *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)&
1096 &                 -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1097                end do
1098              end do
1099 !$OMP END DO
1100            else
1101 !$OMP DO
1102              do ipw=1,npw
1103                ztab(ipw)=czero
1104                do ilmn=1,nlmn
1105                  ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) &
1106 &                 -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1107                end do
1108              end do
1109 !$OMP END DO
1110            end if
1111 
1112 !        ------
1113          else if (choice==5) then ! full derivative w.r.t. k
1114            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1115 !$OMP DO
1116            do ipw=1,npw
1117              ztab(ipw)=czero
1118              do ilmn=1,nlmn
1119                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) &
1120 &               +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1121              end do
1122            end do
1123 !$OMP END DO
1124 
1125 !        ------
1126          else if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi>
1127 !$OMP DO
1128            do ipw=1,npw
1129              ztab(ipw)=czero
1130              do ilmn=1,nlmn
1131                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
1132              end do
1133            end do
1134 !$OMP END DO
1135 
1136 !        ------
1137          else if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi>
1138            ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1139 !$OMP DO
1140            do ipw=1,npw
1141              ztab(ipw)=czero
1142              do ilmn=1,nlmn
1143                ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1144              end do
1145            end do
1146 !$OMP END DO
1147 
1148 !        ------
1149          else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>S<dp/dk_(idir-1)|psi> -
1150 !          <G|dp/dk_(idir-1)>V<dp/dk_(idir+1)|psi>
1151            fdf = ffnl_dir_dat(2*idir-1)
1152            fdb = ffnl_dir_dat(2*idir)
1153 !$OMP DO
1154            do ipw=1,npw
1155              ztab(ipw)=czero
1156              do ilmn=1,nlmn
1157                ztab(ipw)=ztab(ipw) &
1158 &               +ffnl(ipw,fdf,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp) &
1159 &               -ffnl(ipw,fdb,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
1160              end do
1161            end do
1162 !$OMP END DO
1163 
1164 !        ------
1165          else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k
1166 !$OMP DO
1167            do ipw=1,npw
1168              ztab(ipw)=czero
1169              do ilmn=1,nlmn
1170                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp)
1171              end do
1172              ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw)
1173              do ilmn=1,nlmn
1174                ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)
1175              end do
1176            end do
1177 !$OMP END DO
1178 
1179 !        ------
1180          else if (choice==8) then ! full second order derivative w.r.t. k
1181            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1182            ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
1183            ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
1184 !$OMP DO
1185            do ipw=1,npw
1186              ztab(ipw)=czero
1187              do ilmn=1,nlmn
1188                ztab(ipw)=ztab(ipw) &
1189 &               +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx(    gxfacs_(1,  ilmn),    gxfacs_(2,  ilmn),kind=dp)&
1190 &               +ffnl(ipw,1            ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
1191 &               +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)&
1192 &               +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
1193              end do
1194            end do
1195 !$OMP END DO
1196 
1197 !        ------
1198          else if (choice==81) then
1199            ! partial second order derivative w.r.t. k
1200            ! full derivative w.r.t. k1, right derivative w.r.t. k2
1201            !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1202            ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
1203 !$OMP DO
1204            do ipw=1,npw
1205              ztab(ipw)=czero
1206              do ilmn=1,nlmn
1207                ztab(ipw)=ztab(ipw) &
1208 &               +ffnl(ipw,1            ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
1209 &               +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
1210              end do
1211            end do
1212 !$OMP END DO
1213 
1214 !        ------
1215          else
1216 !$OMP WORKSHARE
1217            ztab(:)=czero
1218 !$OMP END WORKSHARE
1219          end if
1220 
1221 
1222 !        ------
1223 !        The OMP WORKSHARE directive doesn't have a good performance with Intel Compiler
1224 !        !$OMP WORKSHARE
1225 !        ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
1226 !        !$OMP END WORKSHARE
1227 !        !$OMP WORKSHARE
1228 !        svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
1229 !        svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:))
1230 !        !$OMP END WORKSHARE
1231 !$OMP DO
1232          do ipw=1,npw
1233            ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp)
1234            svect(1,ipw+ipwshft)=svect(1,ipw+ipwshft)+real(ztab(ipw))
1235            svect(2,ipw+ipwshft)=svect(2,ipw+ipwshft)+aimag(ztab(ipw))
1236          end do
1237 !$OMP END DO
1238 !$OMP END PARALLEL
1239        end if
1240 
1241 !      End loop on atoms
1242      end do
1243 !    End loop on spinors
1244    end do
1245 
1246 !  ==========================================================================
1247  end if
1248 
1249  ABI_DEALLOCATE(ztab)
1250 
1251  if (paw_opt/=3) then
1252    ABI_DEALLOCATE(gxfac_)
1253    if (choice>1) then
1254      ABI_DEALLOCATE(dgxdtfac_)
1255    end if
1256    if (choice==54.or.choice==8.or.choice==81) then
1257      ABI_DEALLOCATE(d2gxdtfac_)
1258    end if
1259  end if
1260  if (paw_opt>=3) then
1261    ABI_DEALLOCATE(gxfacs_)
1262    if (choice>1) then
1263      ABI_DEALLOCATE(dgxdtfacs_)
1264    end if
1265    if (choice==54.or.choice==8.or.choice==81) then
1266      ABI_DEALLOCATE(d2gxdtfacs_)
1267    end if
1268  end if
1269 
1270  DBG_EXIT("COLL")
1271 
1272 #if !defined HAVE_OPENMP
1273 !Fake use of unused variable
1274  if (.false.) write(std_out,*) ipw
1275 #endif
1276 
1277 end subroutine opernlb_ylm