TABLE OF CONTENTS


ABINIT/m_opernla_ylm [ Modules ]

[ Top ] [ Modules ]

NAME

  m_opernla_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_opernla_ylm
27 
28  use defs_basis
29  use defs_abitypes
30  use m_abicore
31  use m_errors
32  use m_xmpi
33 #if defined HAVE_OPENMP
34  use OMP_LIB
35 #endif
36 
37  use m_time,        only : timab
38 
39  implicit none
40 
41  private

ABINIT/opernla_ylm [ Functions ]

[ Top ] [ Functions ]

NAME

 opernla_ylm

FUNCTION

 For a given wave-function |c>, get all projected scalars
 <p_lmn|c> where |p_lmn> are non-local projectors
   With:
   <p_lmn|c>=4pi/sqrt(vol) (i)^l Sum_g[c(g).f_nl(g).Y_lm(g).exp(2pi.i.g.R)]

INPUTS

  choice=chooses possible output:
         if choice>=0: compute projected scalars
         if choice<0: same as choice>0 but use already computed projected scalars
         if ABS(choice)>1, then compute additional quantities:
           2: compute projected scalars and derivatives wrt atm pos.
           3: compute projected scalars and derivatives wrt strains
           23: compute projected scalars, derivatives wrt atm pos. and derivatives wrt strains
           4, 24: compute projected scalars, derivatives wrt atm pos.
                  and 2nd derivatives wrt atm pos.
           5,51,52: compute projected scalars and derivatives wrt wave vector k
           53: compute projected scalars and derivatives wrt wave vector k in direction idir+1 and idir-1
           54: compute projected scalars, deriv. wrt atm pos., deriv. wrt wave vector k
               and 2nd derivatives wrt right wave vector k and atm pos.
           55: compute projected scalars, deriv. strains, deriv. wrt wave vector k
               and 2nd derivatives wrt right wave vector k and strain
           6: compute projected scalars, derivatives wrt atm pos., derivatives wrt strains,
              2nd derivatives wrt 2 strains and derivatives wrt strain and atm pos.
           7: not available
           8: compute projected scalars, derivatives wrt wave vector k
              and 2nd derivatives wrt 2 wave vectors k
  cplex=1 if <p_lmn|c> scalars are real or pure imaginary (equivalent to istwfk>1)
        2 if <p_lmn|c> scalars are complex
  dimffnl=second dimension of ffnl
  ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
  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,signs=2)
                        - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
  indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
  istwf_k=option parameter that describes the storage of wfs
  kpg(npw,nkpg)=(k+G) components          for ikpg=1...3   (if nkpg=3 or 9)
       [(k+G)_a].[(k+G)_b] quantities for ikpg=4...9   (if nkpg=9)
  matblk=dimension of the array ph3d
  mpi_enreg=information about MPI parallelization
  ndgxdt=second dimension of dgxdt
  nd2gxdt=second dimension of d2gxdt
  nincat=number of atoms in the subset here treated
  nkpg=second dimension of array kpg (0, 3 or 9)
  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)
  ph3d(2,npw,matblk)=three-dimensional phase factors
  signs=chooses possible output:
   signs=1: compute derivatives in all directions
   signs=2: compute derivative in direction IDIR only
            compatible only with 1st-order derivatives and "single" derivatives
  ucvol=unit cell volume (bohr^3)
  vect(2,npw*my_nspinor)=starting vector in reciprocal space

OUTPUT

  if (choice>1) dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)=
     gradients of projected scalars wrt coords  (choice=2, 23, 4, 54, 6)
                                    wrt strains (choice=3, 23, 55)
                                    wrt k wave vect. (choice=5, 51, 52, 53, 54, 55, 8)
  if (choice=4, 24, 54, 55, 6, 8) d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)=
     2nd grads of projected scalars wrt 2 coords (choice=4 or 24)
                                    wrt coords & k wave vect. (choice=54)
                                    wrt strains & k wave vect. (choice=55)
                                    wrt coords & strains (choice=6)
                                    wrt 2 strains (choice=6)
                                    wrt 2 k wave vect. (choice=8)
     only compatible with signs=1
  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_d2gxdt(nd2gxdt) = used only when cplex = 1
             cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary

SIDE EFFECTS

  gx(cplex,nlmn,nincat,nspinor)= projected scalars - input if choice<0, output if choice>=0

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

      getcprj,nonlop_ylm

CHILDREN

      timab,xmpi_sum

SOURCE

 146 subroutine opernla_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,d2gxdt,dgxdt,ffnl,gx,&
 147 &       ia3,idir,indlmn,istwf_k,kpg,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpg,nlmn,&
 148 &       nloalg,npw,nspinor,ph3d,signs,ucvol,vect)
 149 
 150 
 151 !This section has been created automatically by the script Abilint (TD).
 152 !Do not modify the following lines by hand.
 153 #undef ABI_FUNC
 154 #define ABI_FUNC 'opernla_ylm'
 155 !End of the abilint section
 156 
 157  implicit none
 158 
 159 !Arguments ------------------------------------
 160 !scalars
 161  integer,intent(in) :: choice,cplex,dimffnl,ia3,idir,istwf_k,matblk,nd2gxdt
 162  integer,intent(in) :: ndgxdt,nincat,nkpg,nlmn,npw,nspinor,signs
 163  real(dp),intent(in) :: ucvol
 164  type(MPI_type),intent(in) :: mpi_enreg
 165 !arrays
 166  integer,intent(in) :: indlmn(6,nlmn),nloalg(3)
 167  integer,intent(out) :: cplex_dgxdt(ndgxdt),cplex_d2gxdt(nd2gxdt)
 168  real(dp),intent(in) :: ffnl(npw,dimffnl,nlmn),kpg(npw,nkpg),ph3d(2,npw,matblk)
 169  real(dp),intent(in) :: vect(:,:)
 170  real(dp),intent(out) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)
 171  real(dp),intent(out) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor),gx(cplex,nlmn,nincat,nspinor)
 172 
 173 !Local variables-------------------------------
 174 !scalars
 175  integer :: choice_,gama,gamb,gamc,gamd,ia,iaph3d
 176  integer :: ierr,il,ilmn,ipw,ipw0,ipwshft,ispinor,jpw,mu,mu0
 177  integer :: mua,mub,ndir,nthreads,nua1,nua2,nub1,nub2
 178  real(dp), parameter :: two_pi2=two_pi*two_pi
 179  real(dp) :: aux_i,aux_i2,aux_i3,aux_i4
 180  real(dp) :: aux_r,aux_r2,aux_r3,aux_r4
 181  real(dp) :: buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5,buffer_i6
 182  real(dp) :: buffer_ia,buffer_ib,buffer_ic,buffer_id,buffer_ie,buffer_if
 183  real(dp) :: buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5,buffer_r6
 184  real(dp) :: buffer_ra,buffer_rb,buffer_rc,buffer_rd,buffer_re,buffer_rf
 185  real(dp) :: kpga,kpgb,kpgc,kpgd,scale,scale2,wt
 186  logical :: check,parity
 187 !arrays
 188  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
 189  integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
 190  integer,parameter :: ffnl_dir_dat(6)=(/2,3,3,1,1,2/)
 191  integer,parameter :: idir1(9)=(/1,1,1,2,2,2,3,3,3/),idir2(9)=(/1,2,3,1,2,3,1,2,3/)
 192  integer :: ffnl_dir(2)
 193  real(dp) :: tsec(2)
 194  real(dp),allocatable :: scali(:),scalr(:),scalari(:,:),scalarr(:,:)
 195 ! *************************************************************************
 196 
 197  if (choice==-1) return
 198 
 199 !Useful variables
 200  choice_=abs(choice)
 201  wt=four_pi/sqrt(ucvol);if (cplex==1) wt=2.d0*wt
 202  ipw0=1;if (istwf_k==2.and.mpi_enreg%me_g0==1) ipw0=2
 203  cplex_dgxdt(:)  = 0 ; if (cplex == 1) cplex_dgxdt(:)  = 1
 204  cplex_d2gxdt(:) = 0 ; if (cplex == 1) cplex_d2gxdt(:) = 1
 205  nthreads=1
 206 
 207 !Check compatibility of options
 208  if (signs==1) then
 209 !  signs=1, almost all choices
 210    check=(choice_==0 .or.choice_==1 .or.choice_==2 .or.choice_==3 .or.choice_==4 .or.&
 211 &   choice_==23.or.choice_==24.or.choice_==5 .or.choice_==51.or.choice_==52.or.&
 212 &   choice_==53.or.choice_==54.or.choice_==55.or.&
 213 &   choice_==6 .or.choice_==8 .or.choice_==81)
 214    ABI_CHECK(check,'BUG: choice not compatible (for signs=1)')
 215  end if
 216  if (signs==2) then
 217 !  signs=2,less choices
 218    check=(choice_== 0.or.choice_== 1.or.choice_== 2.or.choice_==3.or.choice_==5.or.&
 219 &   choice_==23.or.choice_==51.or.choice_==52.or.choice_==54.or.choice_==55.or.&
 220 &   choice_== 8.or.choice_==81)
 221    ABI_CHECK(check,'BUG: signs=2 not compatible with this choice')
 222  end if
 223  if (choice_==3.and.signs==2) then
 224 !  1<=idir<=6 is  required when choice=3 and signs=2
 225    check=(idir>=1.and.idir<=6)
 226    ABI_CHECK(check,'BUG: choice=3 and signs=2 requires 1<=idir<=6')
 227  else if ((choice_==8.or.choice_==81.or.choice_==54).and.signs==2) then
 228 !  1<=idir<=9 is required when choice=8 and signs=2
 229    check=(idir>=1.and.idir<=9)
 230    ABI_CHECK(check,'BUG: choice=8/81 and signs=2 requires 1<=idir<=9')
 231  else
 232 !  signs=2 requires 1<=idir<=3 when choice>1
 233    check=(signs/=2.or.choice<=1.or.(idir>=1.and.idir<=3))
 234    ABI_CHECK(check,'BUG: signs=2 requires 1<=idir<=3')
 235  end if
 236 
 237 #if defined HAVE_OPENMP
 238  nthreads=OMP_GET_NUM_THREADS()
 239 #endif
 240 
 241 !==========================================================================
 242 !========== STANDARD VERSION ==============================================
 243 !==========================================================================
 244 !Parallelization OpenMP on nlmn loops
 245  if (nthreads==1) then
 246 
 247 !  Loop on spinorial components
 248    do ispinor =1,nspinor
 249      ipwshft=(ispinor-1)*npw
 250 
 251 !    Allocate work space
 252      ABI_ALLOCATE(scali,(npw))
 253      ABI_ALLOCATE(scalr,(npw))
 254 
 255 !    Loop on atoms (blocking)
 256      do ia=1,nincat
 257        iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
 258 !      Compute c(g).exp(2pi.i.g.R)
 259        do ipw=ipw0,npw
 260          jpw=ipw+ipwshft
 261          scalr(ipw)=(vect(1,jpw)*ph3d(1,ipw,iaph3d)-vect(2,jpw)*ph3d(2,ipw,iaph3d))
 262          scali(ipw)=(vect(2,jpw)*ph3d(1,ipw,iaph3d)+vect(1,jpw)*ph3d(2,ipw,iaph3d))
 263        end do
 264 
 265        if (ipw0==2) then
 266          scalr(1)=half*vect(1,1+ipwshft)*ph3d(1,1,iaph3d)
 267          scali(1)=half*vect(1,1+ipwshft)*ph3d(2,1,iaph3d)
 268        end if
 269 
 270 !      --------------------------------------------------------------------
 271 !      ALL CHOICES:
 272 !      Accumulate Gx
 273 !      --------------------------------------------------------------------
 274 
 275        if (choice>=0) then
 276          do ilmn=1,nlmn
 277            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 278            scale=wt;if (il>1) scale=-scale
 279            if (cplex==2) then
 280              buffer_r1 = zero ; buffer_i1 = zero
 281              do ipw=1,npw
 282                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)
 283                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)
 284              end do
 285              if (parity) then
 286                gx(1,ilmn,ia,ispinor) = scale*buffer_r1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_i1
 287              else
 288                gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_r1
 289              end if
 290            else
 291              if (parity) then
 292                buffer_r1 =  zero
 293                do ipw=1,npw
 294                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)
 295                end do
 296                gx(1,ilmn,ia,ispinor) = scale*buffer_r1
 297              else
 298                buffer_i1 = zero
 299                do ipw=1,npw
 300                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)
 301                end do
 302                gx(1,ilmn,ia,ispinor) =-scale*buffer_i1
 303              end if
 304            end if
 305          end do
 306        end if
 307 
 308 !      --------------------------------------------------------------------
 309 !      CHOICE= 2, 4, 6, 23, 24, 54  --  SIGNS= 1
 310 !      Accumulate dGxdt --- derivative wrt atm pos. --- for all directions
 311 !      --------------------------------------------------------------------
 312        if ((signs==1).and. &
 313 &       (choice_==2.or.choice_==4.or.choice_==6.or.choice_==23.or.choice_==24.or.choice_==54)) then
 314          mu0=1;if (choice_==23.or.choice_==6) mu0=7
 315          do ilmn=1,nlmn
 316            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 317            scale=wt;if (il>1) scale=-scale
 318            scale2=two_pi*scale
 319            if (cplex==2) then
 320              buffer_r1 = zero ; buffer_r2 = zero
 321              buffer_r3 = zero ; buffer_i1 = zero
 322              buffer_i2 = zero ; buffer_i3 = zero
 323              do ipw=1,npw
 324                aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
 325                aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
 326                buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1)
 327                buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2)
 328                buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3)
 329                buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1)
 330                buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2)
 331                buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3)
 332              end do
 333              if (parity) then
 334                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_r1
 335                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_i1
 336                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
 337                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_i2
 338                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
 339                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_i3
 340              else
 341                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_i1
 342                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_r1
 343                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2
 344                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
 345                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3
 346                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
 347              end if
 348            else
 349              if (parity) then
 350                buffer_r1 = zero ; buffer_r2 = zero ; buffer_r3 = zero
 351                do ipw=1,npw
 352                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
 353                  buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1)
 354                  buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2)
 355                  buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3)
 356                end do
 357                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_r1
 358                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
 359                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
 360              else
 361                buffer_i1 = zero ; buffer_i2 = zero ; buffer_i3 = zero
 362                do ipw=1,npw
 363                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
 364                  buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1)
 365                  buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2)
 366                  buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3)
 367                end do
 368                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_i1
 369                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2
 370                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3
 371              end if
 372            end if
 373          end do
 374        end if
 375 
 376 !      --------------------------------------------------------------------
 377 !      CHOICE= 2  --  SIGNS= 2
 378 !      Accumulate dGxdt --- derivative wrt atm pos. --- for direction IDIR
 379 !      --------------------------------------------------------------------
 380        if ((signs==2).and.(choice_==2)) then
 381          mu0=1
 382          do ilmn=1,nlmn
 383            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 384            scale=wt;if (il>1) scale=-scale
 385            scale2=two_pi*scale
 386            if (cplex==2) then
 387              buffer_r1 = zero ; buffer_i1 = zero
 388              do ipw=1,npw
 389                aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
 390                aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
 391                buffer_r1 = buffer_r1 - aux_i*kpg(ipw,idir)
 392                buffer_i1 = buffer_i1 + aux_r*kpg(ipw,idir)
 393              end do
 394              if (parity) then
 395                dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
 396                dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_i1
 397              else
 398                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
 399                dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
 400              end if
 401            else
 402              if (parity) then
 403                buffer_r1 = zero
 404                do ipw=1,npw
 405                  buffer_r1 = buffer_r1 - scali(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir)
 406                end do
 407                dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
 408              else
 409                buffer_i1 = zero
 410                do ipw=1,npw
 411                  buffer_i1 = buffer_i1 + scalr(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir)
 412                end do
 413                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
 414              end if
 415            end if
 416          end do
 417        end if
 418 
 419 !      --------------------------------------------------------------------
 420 !      CHOICE= 3, 6, 23, 55  -- SIGNS= 1
 421 !      Accumulate dGxdt --- derivative wrt strain --- for all directions
 422 !      --------------------------------------------------------------------
 423        if ((signs==1).and.(choice_==3.or.choice_==6.or.choice_==23.or.choice_==55)) then
 424          mu0=1
 425          do ilmn=1,nlmn
 426            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 427            scale=wt;if (il>1) scale=-scale
 428            scale2=half*scale
 429            if (cplex==2) then
 430              buffer_r1 = zero ; buffer_r2 = zero
 431              buffer_r3 = zero ; buffer_r4 = zero
 432              buffer_r5 = zero ; buffer_r6 = zero
 433              buffer_i1 = zero ; buffer_i2 = zero
 434              buffer_i3 = zero ; buffer_i4 = zero
 435              buffer_i5 = zero ; buffer_i6 = zero
 436              do ipw=1,npw
 437                aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn)
 438                aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn)
 439                aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn)
 440                aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn)
 441                aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn)
 442                aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn)
 443                buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1)
 444                buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2)
 445                buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3)
 446                buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3)
 447                buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3)
 448                buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2)
 449                buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1)
 450                buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2)
 451                buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3)
 452                buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3)
 453                buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3)
 454                buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2)
 455              end do
 456              if (parity) then
 457                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_r1
 458                dgxdt(2,mu0  ,ilmn,ia,ispinor) =-scale*buffer_i1
 459                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
 460                dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2
 461                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
 462                dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3
 463                dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
 464                dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_i4
 465                dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
 466                dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_i5
 467                dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
 468                dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_i6
 469              else
 470                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
 471                dgxdt(2,mu0  ,ilmn,ia,ispinor) =-scale*buffer_r1
 472                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
 473                dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
 474                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
 475                dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
 476                dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4
 477                dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
 478                dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5
 479                dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
 480                dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6
 481                dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
 482              end if
 483            else
 484              if (parity) then
 485                buffer_r1 = zero ; buffer_r2 = zero
 486                buffer_r3 = zero ; buffer_r4 = zero
 487                buffer_r5 = zero ; buffer_r6 = zero
 488                do ipw=1,npw
 489                  aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn)
 490                  aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn)
 491                  aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn)
 492                  buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1)
 493                  buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2)
 494                  buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3)
 495                  buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3)
 496                  buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3)
 497                  buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2)
 498                end do
 499                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_r1
 500                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
 501                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
 502                dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
 503                dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
 504                dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
 505              else
 506                buffer_i1 = zero ; buffer_i2 = zero
 507                buffer_i3 = zero ; buffer_i4 = zero
 508                buffer_i5 = zero ; buffer_i6 = zero
 509                do ipw=1,npw
 510                  aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn)
 511                  aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn)
 512                  aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn)
 513                  buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1)
 514                  buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2)
 515                  buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3)
 516                  buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3)
 517                  buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3)
 518                  buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2)
 519                end do
 520                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
 521                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
 522                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
 523                dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4
 524                dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5
 525                dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6
 526              end if
 527            end if
 528          end do
 529        end if
 530 
 531 !      --------------------------------------------------------------------
 532 !      CHOICE= 3  --  SIGNS= 2
 533 !      Accumulate dGxdt --- derivative wrt strain --- for direction IDIR
 534 !      --------------------------------------------------------------------
 535        if ((signs==2).and.(choice_==3)) then
 536          mu0=1;ffnl_dir(1)=1;if (dimffnl>2) ffnl_dir(1)=idir
 537          do ilmn=1,nlmn
 538            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 539            scale=wt;if (il>1) scale=-scale
 540            if (cplex==2) then
 541              buffer_r1 = zero ; buffer_i1 = zero
 542              do ipw=1,npw
 543                buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
 544                buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
 545              end do
 546              if (parity) then
 547                dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
 548                dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1
 549              else
 550                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
 551                dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1
 552              end if
 553            else
 554              if (parity) then
 555                buffer_r1 = zero
 556                do ipw=1,npw
 557                  buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
 558                end do
 559                dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
 560              else
 561                buffer_i1 = zero
 562                do ipw=1,npw
 563                  buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
 564                end do
 565                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
 566              end if
 567            end if
 568          end do
 569        end if
 570 
 571 !      --------------------------------------------------------------------
 572 !      CHOICE= 4, 24  -- SIGNS= 1
 573 !      Accumulate d2Gxdt --- 2nd derivative wrt 2 atm pos. --- for all dirs
 574 !      --------------------------------------------------------------------
 575        if ((signs==1).and.(choice_==4.or.choice_==24)) then
 576          mu0=1
 577          do ilmn=1,nlmn
 578            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 579            scale=wt;if (il>1) scale=-scale
 580            scale2=two_pi2*scale
 581            if (cplex==2) then
 582              buffer_ra = zero ; buffer_rb = zero
 583              buffer_rc = zero ; buffer_rd = zero
 584              buffer_re = zero ; buffer_rf = zero
 585              buffer_ia = zero ; buffer_ib = zero
 586              buffer_ic = zero ; buffer_id = zero
 587              buffer_ie = zero ; buffer_if = zero
 588              do ipw=1,npw
 589                aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
 590                aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
 591                buffer_ra = buffer_ra - aux_r*kpg(ipw,4)
 592                buffer_rb = buffer_rb - aux_r*kpg(ipw,5)
 593                buffer_rc = buffer_rc - aux_r*kpg(ipw,6)
 594                buffer_rd = buffer_rd - aux_r*kpg(ipw,7)
 595                buffer_re = buffer_re - aux_r*kpg(ipw,8)
 596                buffer_rf = buffer_rf - aux_r*kpg(ipw,9)
 597                buffer_ia = buffer_ia - aux_i*kpg(ipw,4)
 598                buffer_ib = buffer_ib - aux_i*kpg(ipw,5)
 599                buffer_ic = buffer_ic - aux_i*kpg(ipw,6)
 600                buffer_id = buffer_id - aux_i*kpg(ipw,7)
 601                buffer_ie = buffer_ie - aux_i*kpg(ipw,8)
 602                buffer_if = buffer_if - aux_i*kpg(ipw,9)
 603              end do
 604              if (parity) then
 605                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ra
 606                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ia
 607                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb
 608                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_ib
 609                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc
 610                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_ic
 611                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd
 612                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_id
 613                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re
 614                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_ie
 615                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf
 616                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_if
 617              else
 618                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_ia
 619                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ra
 620                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib
 621                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb
 622                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic
 623                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc
 624                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id
 625                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd
 626                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie
 627                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re
 628                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if
 629                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf
 630              end if
 631            else
 632              if (parity) then
 633                buffer_ra = zero ; buffer_rb = zero
 634                buffer_rc = zero ; buffer_rd = zero
 635                buffer_re = zero ; buffer_rf = zero
 636                do ipw=1,npw
 637                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
 638                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
 639                  buffer_ra = buffer_ra - aux_r*kpg(ipw,4)
 640                  buffer_rb = buffer_rb - aux_r*kpg(ipw,5)
 641                  buffer_rc = buffer_rc - aux_r*kpg(ipw,6)
 642                  buffer_rd = buffer_rd - aux_r*kpg(ipw,7)
 643                  buffer_re = buffer_re - aux_r*kpg(ipw,8)
 644                  buffer_rf = buffer_rf - aux_r*kpg(ipw,9)
 645                end do
 646                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ra
 647                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb
 648                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc
 649                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd
 650                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re
 651                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf
 652              else
 653                buffer_ia = zero ; buffer_ib = zero
 654                buffer_ic = zero ; buffer_id = zero
 655                buffer_ie = zero ; buffer_if = zero
 656                do ipw=1,npw
 657                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
 658                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
 659                  buffer_ia = buffer_ia - aux_i*kpg(ipw,4)
 660                  buffer_ib = buffer_ib - aux_i*kpg(ipw,5)
 661                  buffer_ic = buffer_ic - aux_i*kpg(ipw,6)
 662                  buffer_id = buffer_id - aux_i*kpg(ipw,7)
 663                  buffer_ie = buffer_ie - aux_i*kpg(ipw,8)
 664                  buffer_if = buffer_if - aux_i*kpg(ipw,9)
 665                end do
 666                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_ia
 667                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib
 668                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic
 669                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id
 670                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie
 671                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if
 672              end if
 673            end if
 674          end do
 675        end if
 676 
 677 !      --------------------------------------------------------------------
 678 !      CHOICE= 5, 51, 52, 53, 54, 55, 8, 81  --  SIGNS= 1
 679 !      Accumulate dGxdt --- derivative wrt k --- for all directions
 680 !      --------------------------------------------------------------------
 681        if ((signs==1).and.&
 682 &       (choice_==5 .or.choice_==51.or.choice_==52.or.choice_==53.or.&
 683 &       choice_==54.or.choice_==55.or.choice_==8 .or.choice_==81)) then
 684          mu0=1
 685          if (choice_==54) mu0=4
 686          if (choice_==55) mu0=7
 687          do ilmn=1,nlmn
 688            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 689            scale=wt;if (il>1) scale=-scale
 690            if (cplex==2) then
 691              buffer_r1 = zero ; buffer_r2 = zero
 692              buffer_r3 = zero
 693              buffer_i1 = zero ; buffer_i2 = zero
 694              buffer_i3 = zero
 695              do ipw=1,npw
 696                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn)
 697                buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn)
 698                buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn)
 699                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn)
 700                buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn)
 701                buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn)
 702              end do
 703              if (parity) then
 704                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_r1
 705                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
 706                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
 707                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
 708                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
 709                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
 710              else
 711                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_i1
 712                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_r1
 713                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2
 714                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
 715                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3
 716                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
 717              end if
 718            else
 719              cplex_dgxdt(mu0:mu0+2) = 2 ! Warning dgxdt is here pure imaginary
 720              if (parity) then
 721                buffer_i1 = zero ; buffer_i2 = zero
 722                buffer_i3 = zero
 723                do ipw=1,npw
 724                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn)
 725                  buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn)
 726                  buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn)
 727                end do
 728                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
 729                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
 730                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
 731              else
 732                buffer_r1 = zero ; buffer_r2 = zero
 733                buffer_r3 = zero
 734                do ipw=1,npw
 735                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn)
 736                  buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn)
 737                  buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn)
 738                end do
 739                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_r1
 740                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
 741                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
 742              end if
 743            end if
 744          end do
 745        end if
 746 
 747 !      --------------------------------------------------------------------
 748 !      CHOICE= 5, 51, 52, 53, 54, 8, 81  --  SIGNS= 2
 749 !      Accumulate dGxdt --- derivative wrt k ---
 750 !                            for direction IDIR (CHOICE=5,51,52)
 751 !                            for direction IDIR2 (CHOICE=54,81), only if choice > 0
 752 !                            for directions IDIR-1 & IDIR+1 (CHOICE=53)
 753 !                            for 2 directions (CHOICE=8), only if choice > 0
 754 !      --------------------------------------------------------------------
 755        if ((signs==2).and.&
 756 &       (choice_==5.or.choice_==51.or.choice_==52.or.choice_==53.or.choice_==54.or. &
 757 &       choice==8.or.choice==81)) then ! Note the use of choice and not choice_=abs(choice)
 758          mu0=1
 759          if (choice_==5.or.choice_==51.or.choice_==52) then
 760 !          We compute the derivative in one direction
 761            ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir
 762            ndir=1
 763          end if
 764          if (choice_==53) then
 765 !          Case 53 though need multiple directions and here ffnl will contain the
 766 !          derivative information in locations 2,3, and 4 corresponding to idir = 1, 2,
 767 !          and 3. Moreover, choice 53 needs the derivatives in direction idir+1 and idir-1.
 768 !          The parameter vector ffnl_dir_dat contains the necessary translations in
 769 !          locations 1,2 for idir=1; 3,4 for idir=2; and 5,6 for idir=3.
 770            ffnl_dir(1)=ffnl_dir_dat(2*idir-1);ffnl_dir(2)=ffnl_dir_dat(2*idir)
 771            ndir=2
 772          end if
 773          if (choice_==54) then
 774 !          We compute the derivative in one direction (the second one)
 775            ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir2(idir)
 776            ndir=1
 777          end if
 778          if (choice_==8) then
 779 !          We compute the derivative in two directions
 780 !          depending on the input idir argument
 781            ffnl_dir(1)=idir1(idir);ffnl_dir(2)=idir2(idir)
 782            ndir=2
 783          end if
 784          if (choice_==81) then
 785 !          We compute the derivative in one direction (the second one)
 786            ffnl_dir(1)=idir2(idir)
 787            ndir=1
 788          end if
 789          do mua=1,ndir
 790            mu=mu0+mua-1
 791            do ilmn=1,nlmn
 792              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 793              scale=wt;if (il>1) scale=-scale
 794              if (cplex==2) then
 795                buffer_r1 = zero ; buffer_i1 = zero
 796                do ipw=1,npw
 797                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
 798                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
 799                end do
 800                if (parity) then
 801                  dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1
 802                  dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_i1
 803                else
 804                  dgxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1
 805                  dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_r1
 806                end if
 807              else
 808                cplex_dgxdt(mu) = 2 ! Warning dgxdt is here pure imaginary
 809                if (parity) then
 810                  buffer_i1 = zero
 811                  do ipw=1,npw
 812                    buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
 813                  end do
 814                  dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1
 815                else
 816                  buffer_r1 = zero
 817                  do ipw=1,npw
 818                    buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
 819                  end do
 820                  dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1
 821                end if
 822              end if
 823            end do
 824          end do
 825        end if
 826 
 827 !      --------------------------------------------------------------------
 828 !      CHOICE= 54  --  SIGNS= 1
 829 !      Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- for all dirs
 830 !      --------------------------------------------------------------------
 831        if ((signs==1).and.(choice_==54)) then
 832          mu0=1
 833          do ilmn=1,nlmn
 834            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 835            scale=wt;if (il>1) scale=-scale
 836            scale2=two_pi*scale
 837            if (cplex==2) then
 838              do mua=1,3 ! atm. pos
 839                do mub=1,3 ! k
 840                  mu=mu0+(mub-1)+3*(mua-1)
 841                  buffer_r1 = zero ; buffer_i1 = zero
 842                  do ipw=1,npw
 843                    buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 844                    buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 845                  end do
 846                  if (parity) then
 847                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
 848                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1
 849                  else
 850                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_r1
 851                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
 852                  end if
 853                end do
 854              end do
 855            else
 856              cplex_d2gxdt(mu0:mu0+8) = 2 ! Warning d2gxdt is here pure imaginary
 857              if (parity) then
 858                do mua=1,3 ! atm. pos
 859                  do mub=1,3 ! k
 860                    mu=mu0+(mub-1)+3*(mua-1)
 861                    buffer_r1 = zero
 862                    do ipw=1,npw
 863                      buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 864                    end do
 865                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
 866                  end do
 867                end do
 868              else
 869                do mua=1,3 ! atm. pos
 870                  do mub=1,3 ! k
 871                    mu=mu0+(mub-1)+3*(mua-1)
 872                    buffer_i1 = zero
 873                    do ipw=1,npw
 874                      buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 875                    end do
 876                    d2gxdt(1,mu,ilmn,ia,ispinor) = -scale2*buffer_i1
 877                  end do
 878                end do
 879              end if
 880            end if
 881          end do
 882        end if
 883 
 884 !      --------------------------------------------------------------------
 885 !      CHOICE= 54  --  SIGNS= 2
 886 !      Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- for direction IDIR
 887 !      --------------------------------------------------------------------
 888        if ((signs==2).and.(choice_==54)) then
 889          mu0=1
 890          !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
 891          mua=idir1(idir) ! atm. pos
 892          mub=1;if(dimffnl>2) mub=idir2(idir) ! k
 893          do ilmn=1,nlmn
 894            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 895            scale=wt;if (il>1) scale=-scale
 896            scale2=two_pi*scale
 897            if (cplex==2) then
 898              buffer_r1 = zero ; buffer_i1 = zero
 899              do ipw=1,npw
 900                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 901                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 902              end do
 903              if (parity) then
 904                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
 905                d2gxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
 906              else
 907                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_r1
 908                d2gxdt(2,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
 909              end if
 910            else
 911              cplex_d2gxdt(mu0) = 2 ! Warning d2gxdt is here pure imaginary
 912              if (parity) then
 913                buffer_r1 = zero
 914                do ipw=1,npw
 915                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 916                end do
 917                d2gxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
 918              else
 919                buffer_i1 = zero
 920                do ipw=1,npw
 921                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
 922                end do
 923                d2gxdt(1,mu0,ilmn,ia,ispinor) = -scale2*buffer_i1
 924              end if
 925            end if
 926          end do
 927        end if
 928 
 929 !      --------------------------------------------------------------------
 930 !      CHOICE= 55  --  SIGNS= 1
 931 !      Accumulate d2Gxdt --- 2nd derivative wrt k and strains --- for all dirs
 932 !      --------------------------------------------------------------------
 933        if ((signs==1).and.(choice_==55)) then
 934          mu0=1
 935          do ilmn=1,nlmn
 936            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 937            scale=wt;if (il>1) scale=-scale
 938            if (cplex==2) then
 939              do mua=1,6 ! strain
 940                do mub=1,3 ! k
 941                  mu=mu0+(mub-1)+3*(mua-1)
 942                  buffer_r1 = zero ; buffer_i1 = zero
 943                  do ipw=1,npw
 944                    aux_r=ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub)
 945                    buffer_r1 = buffer_r1 + aux_r*scalr(ipw)
 946                    buffer_i1 = buffer_i1 + aux_r*scali(ipw)
 947                  end do
 948                  if (parity) then
 949                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1
 950                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_i1
 951                  else
 952                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1
 953                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_r1
 954                  end if
 955                end do
 956              end do
 957            else
 958              cplex_d2gxdt(mu0:mu0+17) = 2 ! Warning d2gxdt is here pure imaginary
 959              if (parity) then
 960                do mua=1,6 ! strain
 961                  do mub=1,3 ! k
 962                    mu=mu0+(mub-1)+3*(mua-1)
 963                    buffer_i1 = zero
 964                    do ipw=1,npw
 965                      buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+mua,ilmn)*kpg(ipw,mub)
 966                    end do
 967                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1
 968                  end do
 969                end do
 970              else
 971                do mua=1,6 ! strain
 972                  do mub=1,3 ! k
 973                    mu=mu0+(mub-1)+3*(mua-1)
 974                    buffer_r1 = zero
 975                    do ipw=1,npw
 976                      buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+mua,ilmn)*kpg(ipw,mub)
 977                    end do
 978                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1
 979                  end do
 980                end do
 981              end if
 982            end if
 983          end do
 984        end if
 985 
 986 !      --------------------------------------------------------------------
 987 !      CHOICE= 6  --  SIGNS= 1
 988 !      Accumulate d2Gxdt --- 2nd derivative wrt atm. pos and strain --- for all dirs
 989 !        --------------------------------------------------------------------
 990        if ((signs==1).and.(choice_==6)) then
 991          mu0=1;if (choice_==6) mu0=37
 992          ABI_ALLOCATE(scalarr,(npw,4))
 993          ABI_ALLOCATE(scalari,(npw,4))
 994          do ilmn=1,nlmn
 995            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
 996            scale=wt;if (il>1) scale=-scale
 997            scale2=pi*scale
 998            do mu=1,4
 999              do ipw=1,npw
1000                scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn)
1001                scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn)
1002              end do
1003            end do
1004            if(cplex==2) then
1005              do mub=1,6
1006                nub1=alpha(mub);nub2=beta(mub)
1007                do mua=1,3
1008                  mu=mu0+(mua-1)+3*(mub-1)
1009                  buffer_r1 = zero ; buffer_i1 = zero
1010                  do ipw=1,npw
1011                    buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + &
1012 &                   scalarr(ipw,1+nub2)*kpg(ipw,nub1))
1013                    buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + &
1014 &                   scalari(ipw,1+nub2)*kpg(ipw,nub1))
1015                  end do
1016                  if (parity) then
1017                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1
1018                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_r1
1019                  else
1020                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
1021                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1
1022                  end if
1023                end do
1024              end do
1025            else
1026              if (parity) then
1027                do mub=1,6
1028                  nub1=alpha(mub);nub2=beta(mub)
1029                  do mua=1,3
1030                    mu=mu0+(mua-1)+3*(mub-1)
1031                    buffer_i1 = zero
1032                    do ipw=1,npw
1033                      buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + &
1034 &                     scalari(ipw,1+nub2)*kpg(ipw,nub1))
1035                    end do
1036                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1
1037                  end do
1038                end do
1039              else
1040                do mub=1,6
1041                  nub1=alpha(mub);nub2=beta(mub)
1042                  do mua=1,3
1043                    mu=mu0+(mua-1)+3*(mub-1)
1044                    buffer_r1 = zero
1045                    do ipw=1,npw
1046                      buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + &
1047 &                     scalarr(ipw,1+nub2)*kpg(ipw,nub1))
1048                    end do
1049                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
1050                  end do
1051                end do
1052              end if
1053            end if
1054          end do
1055          ABI_DEALLOCATE(scalarr)
1056          ABI_DEALLOCATE(scalari)
1057        end if
1058 
1059 !      --------------------------------------------------------------------
1060 !      CHOICE= 6  --  SIGNS= 1
1061 !      Accumulate d2Gxdt --- 2nd derivative wrt 2 strains --- for all dirs
1062 !        --------------------------------------------------------------------
1063        if ((signs==1).and.(choice_==6)) then
1064          mu0=1
1065          ABI_ALLOCATE(scalarr,(npw,10))
1066          ABI_ALLOCATE(scalari,(npw,10))
1067          do ilmn=1,nlmn
1068            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1069            scale=wt;if (il>1) scale=-scale
1070            scale2=quarter*scale
1071            do mu=1,10
1072              do ipw=1,npw
1073                scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn)
1074                scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn)
1075              end do
1076            end do
1077            if(cplex==2) then
1078              do mub=1,6
1079                nub1=alpha(mub);nub2=beta(mub)
1080                do mua=1,6
1081                  mu=mu0+(mua-1)+6*(mub-1)
1082                  nua1=alpha(mua);nua2=beta(mua)
1083                  gama=gamma(nub1,nua1)
1084                  gamb=gamma(nub2,nua1)
1085                  gamc=gamma(nub1,nua2)
1086                  gamd=gamma(nub2,nua2)
1087                  buffer_r1 = zero ; buffer_i1 = zero
1088                  do ipw=1,npw
1089                    kpga=kpg(ipw,nub1)
1090                    kpgb=kpg(ipw,nub2)
1091                    kpgc=kpg(ipw,nua1)
1092                    kpgd=kpg(ipw,nua2)
1093                    buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd &
1094 &                   + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc
1095                    buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd &
1096 &                   + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc
1097                  end do
1098                  if (parity) then
1099                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
1100                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1
1101                  else
1102                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
1103                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1
1104                  end if
1105                end do
1106              end do
1107            else
1108              if (parity) then
1109                do mub=1,6
1110                  nub1=alpha(mub);nub2=beta(mub)
1111                  do mua=1,6
1112                    mu=mu0+(mua-1)+6*(mub-1)
1113                    nua1=alpha(mua);nua2=beta(mua)
1114                    gama=gamma(nub1,nua1)
1115                    gamb=gamma(nub2,nua1)
1116                    gamc=gamma(nub1,nua2)
1117                    gamd=gamma(nub2,nua2)
1118                    buffer_r1 = zero
1119                    do ipw=1,npw
1120                      kpga=kpg(ipw,nub1)
1121                      kpgb=kpg(ipw,nub2)
1122                      kpgc=kpg(ipw,nua1)
1123                      kpgd=kpg(ipw,nua2)
1124                      buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd &
1125 &                     + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc
1126                    end do
1127                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
1128                  end do
1129                end do
1130              else
1131                do mub=1,6
1132                  nub1=alpha(mub);nub2=beta(mub)
1133                  do mua=1,6
1134                    mu=mu0+(mua-1)+6*(mub-1)
1135                    nua1=alpha(mua);nua2=beta(mua)
1136                    gama=gamma(nub1,nua1)
1137                    gamb=gamma(nub2,nua1)
1138                    gamc=gamma(nub1,nua2)
1139                    gamd=gamma(nub2,nua2)
1140                    buffer_i1 = zero
1141                    do ipw=1,npw
1142                      kpga=kpg(ipw,nub1)
1143                      kpgb=kpg(ipw,nub2)
1144                      kpgc=kpg(ipw,nua1)
1145                      kpgd=kpg(ipw,nua2)
1146                      buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd &
1147 &                     + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc
1148                    end do
1149                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
1150                  end do
1151                end do
1152              end if
1153            end if
1154          end do
1155          ABI_DEALLOCATE(scalarr)
1156          ABI_DEALLOCATE(scalari)
1157        end if
1158 
1159 !      --------------------------------------------------------------------
1160 !      CHOICE= 8, 81  --  SIGNS= 1
1161 !      Accumulate d2Gxdt --- 2nd derivative wrt 2 k wave vect. --- for all dirs
1162 !      --------------------------------------------------------------------
1163        if ((signs==1).and.(choice_==8.or.choice_==81)) then
1164          mu0=1
1165          do ilmn=1,nlmn
1166            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1167            scale=wt;if (il>1) scale=-scale
1168            if (cplex==2) then
1169              buffer_ra = zero ; buffer_rb = zero
1170              buffer_rc = zero ; buffer_rd = zero
1171              buffer_re = zero ; buffer_rf = zero
1172              buffer_ia = zero ; buffer_ib = zero
1173              buffer_ic = zero ; buffer_id = zero
1174              buffer_ie = zero ; buffer_if = zero
1175              do ipw=1,npw
1176                buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn)
1177                buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn)
1178                buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn)
1179                buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn)
1180                buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn)
1181                buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn)
1182                buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn)
1183                buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn)
1184                buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn)
1185                buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn)
1186                buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn)
1187                buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn)
1188              end do
1189              if (parity) then
1190                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_ra
1191                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_ia
1192                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb
1193                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_ib
1194                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc
1195                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_ic
1196                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd
1197                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_id
1198                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re
1199                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_ie
1200                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf
1201                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_if
1202              else
1203                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_ia
1204                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_ra
1205                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib
1206                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb
1207                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic
1208                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc
1209                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id
1210                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd
1211                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie
1212                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_re
1213                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if
1214                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf
1215              end if
1216            else
1217              if (parity) then
1218                buffer_ra = zero ; buffer_rb = zero
1219                buffer_rc = zero ; buffer_rd = zero
1220                buffer_re = zero ; buffer_rf = zero
1221                do ipw=1,npw
1222                  buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn)
1223                  buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn)
1224                  buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn)
1225                  buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn)
1226                  buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn)
1227                  buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn)
1228                end do
1229                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_ra
1230                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb
1231                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc
1232                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd
1233                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re
1234                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf
1235              else
1236                buffer_ia = zero ; buffer_ib = zero
1237                buffer_ic = zero ; buffer_id = zero
1238                buffer_ie = zero ; buffer_if = zero
1239                do ipw=1,npw
1240                  buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn)
1241                  buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn)
1242                  buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn)
1243                  buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn)
1244                  buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn)
1245                  buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn)
1246                end do
1247                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_ia
1248                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib
1249                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic
1250                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id
1251                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie
1252                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if
1253              end if
1254            end if
1255          end do
1256        end if
1257 
1258 !      --------------------------------------------------------------------
1259 !      CHOICE= 8, 81  --  SIGNS= 2
1260 !      Accumulate d2Gxdt --- 2nd derivative wrt k --- for direction IDIR (Voigt not.)
1261 !      --------------------------------------------------------------------
1262        if ((signs==2).and.(choice_==8.or.choice_==81)) then
1263          mu0=1
1264          !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1265          mua=idir1(idir);mub=idir2(idir)
1266          !idir->(Voigt) (11->1,22->2,33->3,32->4,31->5,21->6)
1267          ffnl_dir(1)=gamma(mua,mub)
1268          do ilmn=1,nlmn
1269            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1270            scale=wt;if (il>1) scale=-scale
1271            if (cplex==2) then
1272              buffer_r1 = zero ; buffer_i1 = zero
1273              do ipw=1,npw
1274                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
1275                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
1276              end do
1277              if (parity) then
1278                d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
1279                d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1
1280              else
1281                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
1282                d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1
1283              end if
1284            else
1285              if (parity) then
1286                buffer_r1 = zero
1287                do ipw=1,npw
1288                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
1289                end do
1290                d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
1291              else
1292                buffer_i1 = zero
1293                do ipw=1,npw
1294                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
1295                end do
1296                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
1297              end if
1298            end if
1299          end do
1300        end if
1301 
1302 !      --------------------------------------------------------------------
1303 !      END CHOICES
1304 !      --------------------------------------------------------------------
1305 
1306      end do ! End loop on atoms
1307 
1308 !    Deallocate temporary space
1309      ABI_DEALLOCATE(scali)
1310      ABI_DEALLOCATE(scalr)
1311 
1312    end do !  End loop on spinorial components
1313 
1314 
1315 !  ==========================================================================
1316 !  ========== OPENMP VERSION ================================================
1317 !  ==========================================================================
1318 !  Parallelization OpenMP on npw loops
1319  else
1320 
1321 !$OMP PARALLEL PRIVATE(il,ilmn,ipw,jpw,parity,ndir,ffnl_dir), &
1322 !$OMP PRIVATE(mu,mua,mub,nua1,nua2,nub1,nub2), &
1323 !$OMP PRIVATE(kpga,kpgb,kpgc,kpgd), &
1324 !$OMP PRIVATE(scale,scale2,mu0), &
1325 !$OMP PRIVATE(ispinor,ipwshft,ia,iaph3d)
1326 
1327 !  Loop on spinorial components
1328    do ispinor =1,nspinor
1329      ipwshft=(ispinor-1)*npw
1330 
1331 !    Allocate work space
1332 !$OMP SECTIONS
1333 !$OMP SECTION
1334      ABI_ALLOCATE(scali,(npw))
1335 !$OMP SECTION
1336      ABI_ALLOCATE(scalr,(npw))
1337 !$OMP END SECTIONS
1338 
1339 !    Loop on atoms (blocking)
1340      do ia=1,nincat
1341        iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
1342 
1343 !      Compute Sum_g[c(g).exp(2pi.i.g.R)]
1344 !$OMP DO
1345        do ipw=ipw0,npw
1346          jpw=ipw+ipwshft
1347          scalr(ipw)=(vect(1,jpw)*ph3d(1,ipw,iaph3d)-vect(2,jpw)*ph3d(2,ipw,iaph3d))
1348          scali(ipw)=(vect(2,jpw)*ph3d(1,ipw,iaph3d)+vect(1,jpw)*ph3d(2,ipw,iaph3d))
1349        end do
1350 !$OMP END DO
1351 !$OMP SINGLE
1352        if (ipw0==2) then
1353          scalr(1)=half*vect(1,1+ipwshft)*ph3d(1,1,iaph3d)
1354          scali(1)=half*vect(1,1+ipwshft)*ph3d(2,1,iaph3d)
1355        end if
1356 !$OMP END SINGLE
1357 
1358 !      --------------------------------------------------------------------
1359 !      ALL CHOICES
1360 !      Accumulate Gx
1361 !      --------------------------------------------------------------------
1362        if (choice>=0) then ! JWZ to check: I dont think 53 needs this
1363          do ilmn=1,nlmn
1364            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1365            scale=wt;if (il>1) scale=-scale
1366            if (cplex==2) then
1367 !$OMP SINGLE
1368              buffer_r1 = zero ; buffer_i1 = zero
1369 !$OMP END SINGLE
1370 !$OMP DO &
1371 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
1372              do ipw=1,npw
1373                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)
1374                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)
1375              end do
1376 !$OMP END DO
1377 !$OMP SINGLE
1378              if (parity) then
1379                gx(1,ilmn,ia,ispinor) = scale*buffer_r1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_i1
1380              else
1381                gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_r1
1382              end if
1383 !$OMP END SINGLE
1384            else
1385              if (parity) then
1386 !$OMP SINGLE
1387                buffer_r1 =  zero
1388 !$OMP END SINGLE
1389 !$OMP DO &
1390 !$OMP REDUCTION(+:buffer_r1)
1391                do ipw=1,npw
1392                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)
1393                end do
1394 !$OMP END DO
1395 !$OMP SINGLE
1396                gx(1,ilmn,ia,ispinor) = scale*buffer_r1
1397 !$OMP END SINGLE
1398              else
1399 !$OMP SINGLE
1400                buffer_i1 = zero
1401 !$OMP END SINGLE
1402 !$OMP DO &
1403 !$OMP REDUCTION(+:buffer_i1)
1404                do ipw=1,npw
1405                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)
1406                end do
1407 !$OMP END DO
1408 !$OMP SINGLE
1409                gx(1,ilmn,ia,ispinor) =-scale*buffer_i1
1410 !$OMP END SINGLE
1411              end if
1412            end if
1413          end do
1414        end if
1415 
1416 !      --------------------------------------------------------------------
1417 !      CHOICE= 2, 4, 6, 23, 24, 54  --  SIGNS= 1
1418 !      Accumulate dGxdt --- derivative wrt atm pos. --- for all directions
1419 !      --------------------------------------------------------------------
1420        if ((signs==1).and. &
1421 &       (choice_==2.or.choice_==23.or.choice_==24.or.choice_==4.or.choice_==54.or.choice_==6)) then
1422          mu0 =1;if (choice_==23.or.choice_==6) mu0=7
1423          do ilmn=1,nlmn
1424            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1425            scale=wt;if (il>1) scale=-scale
1426            scale2=two_pi*scale
1427            if (cplex==2) then
1428 !$OMP SINGLE
1429              buffer_r1 = zero ; buffer_r2 = zero
1430              buffer_r3 = zero ; buffer_i1 = zero
1431              buffer_i2 = zero ; buffer_i3 = zero
1432 !$OMP END SINGLE
1433 !$OMP DO &
1434 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_i1,buffer_i2,buffer_i3),&
1435 !$OMP PRIVATE(aux_r,aux_i)
1436              do ipw=1,npw
1437                aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1438                aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1439                buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1)
1440                buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2)
1441                buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3)
1442                buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1)
1443                buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2)
1444                buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3)
1445              end do
1446 !$OMP END DO
1447 !$OMP SINGLE
1448              if (parity) then
1449                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_r1
1450                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_i1
1451                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
1452                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_i2
1453                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
1454                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_i3
1455              else
1456                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_i1
1457                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_r1
1458                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2
1459                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
1460                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3
1461                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
1462              end if
1463 !$OMP END SINGLE
1464            else
1465              if (parity) then
1466 !$OMP SINGLE
1467                buffer_r1 = zero ; buffer_r2 = zero ; buffer_r3 = zero
1468 !$OMP END SINGLE
1469 !$OMP DO &
1470 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3), &
1471 !$OMP PRIVATE(aux_r,aux_i)
1472                do ipw=1,npw
1473                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1474                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1475                  buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1)
1476                  buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2)
1477                  buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3)
1478                end do
1479 !$OMP END DO
1480 !$OMP SINGLE
1481                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_r1
1482                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
1483                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
1484 !$OMP END SINGLE
1485              else
1486 !$OMP SINGLE
1487                buffer_i1 = zero ; buffer_i2 = zero ; buffer_i3 = zero
1488 !$OMP END SINGLE
1489 !$OMP DO &
1490 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3), &
1491 !$OMP PRIVATE(aux_r,aux_i)
1492                do ipw=1,npw
1493                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1494                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1495                  buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1)
1496                  buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2)
1497                  buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3)
1498                end do
1499 !$OMP SINGLE
1500                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_i1
1501                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2
1502                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3
1503 !$OMP END SINGLE
1504              end if
1505            end if
1506          end do
1507        end if
1508 
1509 !      --------------------------------------------------------------------
1510 !      CHOICE= 2  --  SIGNS= 2
1511 !      Accumulate dGxdt --- derivative wrt atm pos. --- for direction IDIR
1512 !      --------------------------------------------------------------------
1513        if ((signs==2).and.(choice_==2)) then
1514          mu0=1
1515          do ilmn=1,nlmn
1516            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1517            scale=wt;if (il>1) scale=-scale
1518            scale2=two_pi*scale
1519            if (cplex==2) then
1520 !$OMP SINGLE
1521              buffer_r1 = zero ; buffer_i1 = zero
1522 !$OMP END SINGLE
1523 !$OMP DO &
1524 !$OMP REDUCTION(+:buffer_r1,buffer_i1), &
1525 !$OMP PRIVATE(aux_r,aux_i)
1526              do ipw=1,npw
1527                aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1528                aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1529                buffer_r1 = buffer_r1 - aux_i*kpg(ipw,idir)
1530                buffer_i1 = buffer_i1 + aux_r*kpg(ipw,idir)
1531              end do
1532 !$OMP END DO
1533 !$OMP SINGLE
1534              if (parity) then
1535                dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
1536                dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_i1
1537              else
1538                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
1539                dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
1540              end if
1541 !$OMP END SINGLE
1542            else
1543              if (parity) then
1544 !$OMP SINGLE
1545                buffer_r1 = zero
1546 !$OMP END SINGLE
1547 !$OMP DO &
1548 !$OMP REDUCTION(+:buffer_r1)
1549                do ipw=1,npw
1550                  buffer_r1 = buffer_r1 - scali(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir)
1551                end do
1552 !$OMP END DO
1553 !$OMP SINGLE
1554                dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
1555 !$OMP END SINGLE
1556              else
1557 !$OMP SINGLE
1558                buffer_i1 = zero
1559 !$OMP END SINGLE
1560 !$OMP DO &
1561 !$OMP REDUCTION(+:buffer_i1)
1562                do ipw=1,npw
1563                  buffer_i1 = buffer_i1 + scalr(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir)
1564                end do
1565 !$OMP END DO
1566 !$OMP SINGLE
1567                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
1568 !$OMP END SINGLE
1569              end if
1570            end if
1571          end do
1572        end if
1573 
1574 !      --------------------------------------------------------------------
1575 !      CHOICE= 3, 6, 23, 55  -- SIGNS= 1
1576 !      Accumulate dGxdt --- derivative wrt strain --- for all directions
1577 !      --------------------------------------------------------------------
1578        if ((signs==1).and.(choice_==3.or.choice_==6.or.choice_==23.or.choice_==55)) then
1579          mu0=1
1580          do ilmn=1,nlmn
1581            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1582            scale=wt;if (il>1) scale=-scale
1583            scale2=half*scale
1584            if (cplex==2) then
1585 !$OMP SINGLE
1586              buffer_r1 = zero ; buffer_r2 = zero
1587              buffer_r3 = zero ; buffer_r4 = zero
1588              buffer_r5 = zero ; buffer_r6 = zero
1589              buffer_i1 = zero ; buffer_i2 = zero
1590              buffer_i3 = zero ; buffer_i4 = zero
1591              buffer_i5 = zero ; buffer_i6 = zero
1592 !$OMP END SINGLE
1593 !$OMP DO &
1594 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5,buffer_r6), &
1595 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5,buffer_i6), &
1596 !$OMP PRIVATE(aux_r2,aux_r3,aux_r4,aux_i2,aux_i3,aux_i4)
1597              do ipw=1,npw
1598                aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn)
1599                aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn)
1600                aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn)
1601                aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn)
1602                aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn)
1603                aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn)
1604                buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1)
1605                buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2)
1606                buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3)
1607                buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3)
1608                buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3)
1609                buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2)
1610                buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1)
1611                buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2)
1612                buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3)
1613                buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3)
1614                buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3)
1615                buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2)
1616              end do
1617 !$OMP END DO
1618 !$OMP SINGLE
1619              if (parity) then
1620                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_r1
1621                dgxdt(2,mu0  ,ilmn,ia,ispinor) =-scale*buffer_i1
1622                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
1623                dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2
1624                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
1625                dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3
1626                dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
1627                dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_i4
1628                dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
1629                dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_i5
1630                dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
1631                dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_i6
1632              else
1633                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
1634                dgxdt(2,mu0  ,ilmn,ia,ispinor) =-scale*buffer_r1
1635                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
1636                dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
1637                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
1638                dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
1639                dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4
1640                dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
1641                dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5
1642                dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
1643                dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6
1644                dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
1645              end if
1646 !$OMP END SINGLE
1647            else
1648              if (parity) then
1649 !$OMP SINGLE
1650                buffer_r1 = zero ; buffer_r2 = zero
1651                buffer_r3 = zero ; buffer_r4 = zero
1652                buffer_r5 = zero ; buffer_r6 = zero
1653 !$OMP END SINGLE
1654 !$OMP DO &
1655 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5,buffer_r6), &
1656 !$OMP PRIVATE(aux_r2,aux_r3,aux_r4)
1657                do ipw=1,npw
1658                  aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn)
1659                  aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn)
1660                  aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn)
1661                  buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1)
1662                  buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2)
1663                  buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3)
1664                  buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3)
1665                  buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3)
1666                  buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2)
1667                end do
1668 !$OMP SINGLE
1669                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_r1
1670                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
1671                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
1672                dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
1673                dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
1674                dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
1675 !$OMP END SINGLE
1676              else
1677 !$OMP SINGLE
1678                buffer_i1 = zero ; buffer_i2 = zero
1679                buffer_i3 = zero ; buffer_i4 = zero
1680                buffer_i5 = zero ; buffer_i6 = zero
1681 !$OMP END SINGLE
1682 !$OMP DO &
1683 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5,buffer_i6), &
1684 !$OMP PRIVATE(aux_i2,aux_i3,aux_i4)
1685                do ipw=1,npw
1686                  aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn)
1687                  aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn)
1688                  aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn)
1689                  buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1)
1690                  buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2)
1691                  buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3)
1692                  buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3)
1693                  buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3)
1694                  buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2)
1695                end do
1696 !$OMP END DO
1697 !$OMP SINGLE
1698                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
1699                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
1700                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
1701                dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4
1702                dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5
1703                dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6
1704 !$OMP END SINGLE
1705              end if
1706            end if
1707          end do
1708        end if
1709 
1710 !      --------------------------------------------------------------------
1711 !      CHOICE= 3  --  SIGNS= 2
1712 !      Accumulate dGxdt --- derivative wrt strain --- for direction IDIR
1713 !      --------------------------------------------------------------------
1714        if ((signs==2).and.(choice_==3)) then
1715          mu0=1;ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir
1716          do ilmn=1,nlmn
1717            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1718            scale=wt;if (il>1) scale=-scale
1719            if (cplex==2) then
1720 !$OMP SINGLE
1721              buffer_r1 = zero ; buffer_i1 = zero
1722 !$OMP END SINGLE
1723 !$OMP DO &
1724 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
1725              do ipw=1,npw
1726                buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
1727                buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
1728              end do
1729 !$OMP END DO
1730 !$OMP SINGLE
1731              if (parity) then
1732                dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
1733                dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1
1734              else
1735                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
1736                dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1
1737              end if
1738 !$OMP END SINGLE
1739            else
1740              if (parity) then
1741 !$OMP SINGLE
1742                buffer_r1 = zero
1743 !$OMP END SINGLE
1744 !$OMP DO &
1745 !$OMP REDUCTION(+:buffer_r1)
1746                do ipw=1,npw
1747                  buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
1748                end do
1749 !$OMP END DO
1750 !$OMP SINGLE
1751                dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
1752 !$OMP END SINGLE
1753              else
1754 !$OMP SINGLE
1755                buffer_i1 = zero
1756 !$OMP END SINGLE
1757 !$OMP DO &
1758 !$OMP REDUCTION(+:buffer_i1)
1759                do ipw=1,npw
1760                  buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn)
1761                end do
1762 !$OMP END DO
1763 !$OMP SINGLE
1764                dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
1765 !$OMP END SINGLE
1766              end if
1767            end if
1768          end do
1769        end if
1770 
1771 !      --------------------------------------------------------------------
1772 !      CHOICE= 4, 24  --  SIGNS= 1
1773 !      Accumulate d2Gxdt --- 2nd derivative wrt 2 atm pos. --- for all dirs
1774 !      --------------------------------------------------------------------
1775        if ((signs==1).and.(choice_==4.or.choice_==24)) then
1776          mu0=1
1777          do ilmn=1,nlmn
1778            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1779            scale=wt;if (il>1) scale=-scale
1780            scale2=two_pi2*scale
1781            if (cplex==2) then
1782 !$OMP SINGLE
1783              buffer_ra = zero ; buffer_rb = zero
1784              buffer_rc = zero ; buffer_rd = zero
1785              buffer_re = zero ; buffer_rf = zero
1786              buffer_ia = zero ; buffer_ib = zero
1787              buffer_ic = zero ; buffer_id = zero
1788              buffer_ie = zero ; buffer_if = zero
1789 !$OMP END SINGLE
1790 !$OMP DO &
1791 !$OMP PRIVATE(aux_r,aux_i), &
1792 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc,buffer_rd,buffer_re,buffer_rf),&
1793 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic,buffer_id,buffer_ie,buffer_if)
1794              do ipw=1,npw
1795                aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1796                aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1797                buffer_ra = buffer_ra - aux_r*kpg(ipw,4)
1798                buffer_rb = buffer_rb - aux_r*kpg(ipw,5)
1799                buffer_rc = buffer_rc - aux_r*kpg(ipw,6)
1800                buffer_rd = buffer_rd - aux_r*kpg(ipw,7)
1801                buffer_re = buffer_re - aux_r*kpg(ipw,8)
1802                buffer_rf = buffer_rf - aux_r*kpg(ipw,9)
1803                buffer_ia = buffer_ia - aux_i*kpg(ipw,4)
1804                buffer_ib = buffer_ib - aux_i*kpg(ipw,5)
1805                buffer_ic = buffer_ic - aux_i*kpg(ipw,6)
1806                buffer_id = buffer_id - aux_i*kpg(ipw,7)
1807                buffer_ie = buffer_ie - aux_i*kpg(ipw,8)
1808                buffer_if = buffer_if - aux_i*kpg(ipw,9)
1809              end do
1810 !$OMP END DO
1811 !$OMP SINGLE
1812              if (parity) then
1813                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ra
1814                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ia
1815                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb
1816                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_ib
1817                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc
1818                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_ic
1819                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd
1820                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_id
1821                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re
1822                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_ie
1823                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf
1824                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_if
1825              else
1826                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_ia
1827                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ra
1828                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib
1829                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb
1830                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic
1831                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc
1832                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id
1833                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd
1834                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie
1835                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re
1836                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if
1837                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf
1838              end if
1839 !$OMP END SINGLE
1840            else
1841              if (parity) then
1842 !$OMP SINGLE
1843                buffer_ra = zero ; buffer_rb = zero
1844                buffer_rc = zero ; buffer_rd = zero
1845                buffer_re = zero ; buffer_rf = zero
1846 !$OMP END SINGLE
1847 !$OMP DO &
1848 !$OMP PRIVATE(aux_r,aux_i), &
1849 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc,buffer_rd,buffer_re,buffer_rf)
1850                do ipw=1,npw
1851                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1852                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1853                  buffer_ra = buffer_ra - aux_r*kpg(ipw,4)
1854                  buffer_rb = buffer_rb - aux_r*kpg(ipw,5)
1855                  buffer_rc = buffer_rc - aux_r*kpg(ipw,6)
1856                  buffer_rd = buffer_rd - aux_r*kpg(ipw,7)
1857                  buffer_re = buffer_re - aux_r*kpg(ipw,8)
1858                  buffer_rf = buffer_rf - aux_r*kpg(ipw,9)
1859                end do
1860 !$OMP END DO
1861 !$OMP SINGLE
1862                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale2*buffer_ra
1863                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb
1864                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc
1865                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd
1866                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re
1867                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf
1868 !$OMP END SINGLE
1869              else
1870 !$OMP SINGLE
1871                buffer_ia = zero ; buffer_ib = zero
1872                buffer_ic = zero ; buffer_id = zero
1873                buffer_ie = zero ; buffer_if = zero
1874 !$OMP END SINGLE
1875 !$OMP DO &
1876 !$OMP PRIVATE(aux_i,aux_r), &
1877 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic,buffer_id,buffer_ie,buffer_if)
1878                do ipw=1,npw
1879                  aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
1880                  aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
1881                  buffer_ia = buffer_ia - aux_i*kpg(ipw,4)
1882                  buffer_ib = buffer_ib - aux_i*kpg(ipw,5)
1883                  buffer_ic = buffer_ic - aux_i*kpg(ipw,6)
1884                  buffer_id = buffer_id - aux_i*kpg(ipw,7)
1885                  buffer_ie = buffer_ie - aux_i*kpg(ipw,8)
1886                  buffer_if = buffer_if - aux_i*kpg(ipw,9)
1887                end do
1888 !$OMP END DO
1889 !$OMP SINGLE
1890                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale2*buffer_ia
1891                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib
1892                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic
1893                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id
1894                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie
1895                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if
1896 !$OMP END SINGLE
1897              end if
1898            end if
1899          end do
1900        end if
1901 
1902 !      --------------------------------------------------------------------
1903 !      CHOICE= 5, 51, 52, 53, 54, 55, 8, 81  --  SIGNS= 1
1904 !      Accumulate dGxdt --- derivative wrt k --- for all directions
1905 !      --------------------------------------------------------------------
1906        if ((signs==1).and.&
1907 &       (choice_==5 .or.choice_==51.or.choice_==52.or.choice_==53.or.&
1908 &       choice_==54.or.choice_==55.or.choice_==8 .or.choice_==81)) then
1909          mu0=1
1910          if (choice_==54) mu0=4
1911          if (choice_==55) mu0=7
1912          do ilmn=1,nlmn
1913            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
1914            scale=wt;if (il>1) scale=-scale
1915            if (cplex==2) then
1916 !$OMP SINGLE
1917              buffer_r1 = zero ; buffer_r2 = zero
1918              buffer_r3 = zero
1919              buffer_i1 = zero ; buffer_i2 = zero
1920              buffer_i3 = zero
1921 !$OMP END SINGLE
1922 !$OMP DO &
1923 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3), &
1924 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3)
1925              do ipw=1,npw
1926                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn)
1927                buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn)
1928                buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn)
1929                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn)
1930                buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn)
1931                buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn)
1932              end do
1933 !$OMP END DO
1934 !$OMP SINGLE
1935              if (parity) then
1936                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_r1
1937                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
1938                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
1939                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
1940                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
1941                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
1942              else
1943                dgxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_i1
1944                dgxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_r1
1945                dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2
1946                dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
1947                dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3
1948                dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
1949              end if
1950 !$OMP END SINGLE
1951            else
1952              cplex_dgxdt(mu0:mu0+2) = 2 ! Warning dgxdt is here pure imaginary
1953              if (parity) then
1954 !$OMP SINGLE
1955                buffer_i1 = zero ; buffer_i2 = zero
1956                buffer_i3 = zero
1957 !$OMP END SINGLE
1958 !$OMP DO &
1959 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3)
1960                do ipw=1,npw
1961                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn)
1962                  buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn)
1963                  buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn)
1964                end do
1965 !$OMP END DO
1966 !$OMP SINGLE
1967                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_i1
1968                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
1969                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
1970 !$OMP END SINGLE
1971              else
1972 !$OMP SINGLE
1973                buffer_r1 = zero ; buffer_r2 = zero
1974                buffer_r3 = zero
1975 !$OMP END SINGLE
1976 !$OMP DO &
1977 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3)
1978                do ipw=1,npw
1979                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn)
1980                  buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn)
1981                  buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn)
1982                end do
1983 !$OMP END DO
1984 !$OMP SINGLE
1985                dgxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_r1
1986                dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
1987                dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
1988 !$OMP END SINGLE
1989              end if
1990            end if
1991          end do
1992        end if
1993 
1994 !      --------------------------------------------------------------------
1995 !      CHOICE= 5, 51, 52, 53, 54, 8, 81  --  SIGNS= 2
1996 !      Accumulate dGxdt --- derivative wrt k ---
1997 !                            for direction IDIR (CHOICE=5,51,52)
1998 !                            for direction IDIR2 (CHOICE=54,81)
1999 !                            for directions IDIR-1 & IDIR+1 (CHOICE=53)
2000 !                            for 2 directions (CHOICE=8)
2001 !      --------------------------------------------------------------------
2002        if ((signs==2).and.&
2003 &       (choice_==5.or.choice_==51.or.choice_==52.or.choice_==53.or.choice_==54.or. &
2004 &       choice_==8.or.choice_==81)) then
2005          mu0=1
2006          if (choice_==5.or.choice_==51.or.choice_==52) then
2007 !          We compute the derivative in one direction
2008            ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir
2009            ndir=1
2010          end if
2011          if (choice_==53) then
2012 !          Case 53 though need multiple directions and here ffnl will contain the
2013 !          derivative information in locations 2,3, and 4 corresponding to idir = 1, 2,
2014 !          and 3. Moreover, choice 53 needs the derivatives in direction idir+1 and idir-1.
2015 !          The parameter vector ffnl_dir_dat contains the necessary translations in
2016 !          locations 1,2 for idir=1; 3,4 for idir=2; and 5,6 for idir=3.
2017            ffnl_dir(1)=ffnl_dir_dat(2*idir-1);ffnl_dir(2)=ffnl_dir_dat(2*idir)
2018            ndir=2
2019          end if
2020          if (choice_==54) then
2021 !          We compute the derivative in one direction (the second one)
2022            ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir2(idir)
2023            ndir=1
2024          end if
2025          if (choice_==8) then
2026 !          We compute the derivative in two directions
2027 !          depending on the input idir argument
2028            ffnl_dir(1)=idir1(idir);ffnl_dir(2)=idir2(idir)
2029            ndir=2
2030          end if
2031          if (choice_==81) then
2032 !          We compute the derivative in one direction (the second one)
2033            ffnl_dir(1)=idir2(idir)
2034            ndir=1
2035          end if
2036          do mua=1,ndir
2037            mu=mu0+mua-1
2038            do ilmn=1,nlmn
2039              il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2040              scale=wt;if (il>1) scale=-scale
2041              if (cplex==2) then
2042 !$OMP SINGLE
2043                buffer_r1 = zero ; buffer_i1 = zero
2044 !$OMP END SINGLE
2045 !$OMP DO &
2046 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
2047                do ipw=1,npw
2048                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
2049                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
2050                end do
2051 !$OMP END DO
2052 !$OMP SINGLE
2053                if (parity) then
2054                  dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1
2055                  dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_i1
2056                else
2057                  dgxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1
2058                  dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_r1
2059                end if
2060 !$OMP END SINGLE
2061              else
2062 !$OMP SINGLE
2063                cplex_dgxdt(mu) = 2 ! Warning dgxdt is here pure imaginary
2064 !$OMP END SINGLE
2065                if (parity) then
2066 !$OMP SINGLE
2067                  buffer_i1 = zero
2068 !$OMP END SINGLE
2069 !$OMP DO &
2070 !$OMP REDUCTION(+:buffer_i1)
2071                  do ipw=1,npw
2072                    buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
2073                  end do
2074 !$OMP END DO
2075 !$OMP SINGLE
2076                  dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1
2077 !$OMP END SINGLE
2078                else
2079 !$OMP SINGLE
2080                  buffer_r1 = zero
2081 !$OMP END SINGLE
2082 !$OMP DO &
2083 !$OMP REDUCTION(+:buffer_r1)
2084                  do ipw=1,npw
2085                    buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
2086                  end do
2087 !$OMP SINGLE
2088                  dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1
2089 !$OMP END SINGLE
2090                end if
2091              end if
2092            end do
2093          end do
2094        end if
2095 
2096 !      --------------------------------------------------------------------
2097 !      CHOICE= 54  --  SIGNS= 1
2098 !      Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- in all dirs
2099 !      --------------------------------------------------------------------
2100        if ((signs==1).and.(choice_==54)) then
2101          mu0=1
2102          do ilmn=1,nlmn
2103            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2104            scale=wt;if (il>1) scale=-scale
2105            scale2=two_pi*scale
2106            if (cplex==2) then
2107              do mua=1,3 ! atm. pos
2108                do mub=1,3 ! k
2109                  mu=mu0+(mub-1)+3*(mua-1)
2110 !$OMP SINGLE
2111                  buffer_r1 = zero ; buffer_i1 = zero
2112 !$OMP END SINGLE
2113 !$OMP DO &
2114 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
2115                  do ipw=1,npw
2116                    buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2117                    buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2118                  end do
2119 !$OMP END DO
2120 !$OMP SINGLE
2121                  if (parity) then
2122                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
2123                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2124                  else
2125                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_r1
2126                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
2127                  end if
2128 !$OMP END SINGLE
2129                end do
2130              end do
2131            else
2132 !$OMP SINGLE
2133              cplex_d2gxdt(mu0:mu0+8) = 2 ! Warning dgxdt is here pure imaginary
2134 !$OMP END SINGLE
2135              if (parity) then
2136                do mua=1,3 ! atm. pos
2137                  do mub=1,3 ! k
2138                    mu=mu0+(mub-1)+3*(mua-1)
2139 !$OMP SINGLE
2140                    buffer_r1 = zero
2141 !$OMP END SINGLE
2142 !$OMP DO &
2143 !$OMP REDUCTION(+:buffer_r1)
2144                    do ipw=1,npw
2145                      buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2146                    end do
2147 !$OMP END DO
2148 !$OMP SINGLE
2149                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2150 !$OMP END SINGLE
2151                  end do
2152                end do
2153              else
2154                do mua=1,3 ! atm. pos
2155                  do mub=1,3 ! k
2156                    mu=mu0+(mub-1)+3*(mua-1)
2157 !$OMP SINGLE
2158                    buffer_i1 = zero
2159 !$OMP END SINGLE
2160 !$OMP DO &
2161 !$OMP REDUCTION(+:buffer_i1)
2162                    do ipw=1,npw
2163                      buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2164                    end do
2165 !$OMP END DO
2166 !$OMP SINGLE
2167                    d2gxdt(1,mu,ilmn,ia,ispinor) = -scale2*buffer_i1
2168 !$OMP END SINGLE
2169                  end do
2170                end do
2171              end if
2172            end if
2173          end do
2174        end if
2175 
2176 !      --------------------------------------------------------------------
2177 !      CHOICE= 54  --  SIGNS= 2
2178 !      Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- for direction IDIR
2179 !      --------------------------------------------------------------------
2180        if ((signs==2).and.(choice_==54)) then
2181          mu0=1
2182          !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
2183          mua=idir1(idir) ! atm. pos
2184          mub=1;if(dimffnl>2) mub=idir2(idir) ! k
2185          do ilmn=1,nlmn
2186            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2187            scale=wt;if (il>1) scale=-scale
2188            scale2=two_pi*scale
2189            if (cplex==2) then
2190 !$OMP SINGLE
2191              buffer_r1 = zero ; buffer_i1 = zero
2192 !$OMP END SINGLE
2193 !$OMP DO &
2194 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
2195              do ipw=1,npw
2196                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2197                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2198              end do
2199 !$OMP END DO
2200 !$OMP SINGLE
2201              if (parity) then
2202                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
2203                d2gxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1
2204              else
2205                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_r1
2206                d2gxdt(2,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1
2207              end if
2208 !$OMP END SINGLE
2209            else
2210              cplex_d2gxdt(mu0) = 2 ! Warning d2gxdt is here pure imaginary
2211              if (parity) then
2212 !$OMP SINGLE
2213                buffer_r1 = zero
2214 !$OMP END SINGLE
2215 !$OMP DO &
2216 !$OMP REDUCTION(+:buffer_r1)
2217                do ipw=1,npw
2218                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2219                end do
2220 !$OMP END DO
2221 !$OMP SINGLE
2222                d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
2223 !$OMP END SINGLE
2224              else
2225 !$OMP SINGLE
2226                buffer_i1 = zero
2227 !$OMP END SINGLE
2228 !$OMP DO &
2229 !$OMP REDUCTION(+:buffer_i1)
2230                do ipw=1,npw
2231                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua)
2232                end do
2233 !$OMP SINGLE
2234                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
2235 !$OMP END SINGLE
2236              end if
2237            end if
2238          end do
2239        end if
2240 
2241 !      --------------------------------------------------------------------
2242 !      CHOICE= 55  --  SIGNS= 1
2243 !      Accumulate d2Gxdt --- 2nd derivative wrt k and strains --- for all dirs
2244 !      --------------------------------------------------------------------
2245        if ((signs==1).and.(choice_==55)) then
2246          mu0=1
2247          do ilmn=1,nlmn
2248            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2249            scale=wt;if (il>1) scale=-scale
2250            if (cplex==2) then
2251              do mua=1,6 ! strain
2252                do mub=1,3 ! k
2253                  mu=mu0+(mub-1)+3*(mua-1)
2254 !$OMP SINGLE
2255                  buffer_r1 = zero ; buffer_i1 = zero
2256 !$OMP END SINGLE
2257 !$OMP DO &
2258 !$OMP REDUCTION(+:buffer_r1,buffer_i1),&
2259 !$OMP PRIVATE(aux_r)
2260                  do ipw=1,npw
2261                    aux_r=ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub)
2262                    buffer_r1 = buffer_r1 + aux_r*scalr(ipw)
2263                    buffer_i1 = buffer_i1 + aux_r*scali(ipw)
2264                  end do
2265 !$OMP END DO
2266 !$OMP SINGLE
2267                  if (parity) then
2268                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1
2269                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_i1
2270                  else
2271                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1
2272                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_r1
2273                  end if
2274 !$OMP END SINGLE
2275                end do
2276              end do
2277            else
2278 !$OMP SINGLE
2279              cplex_d2gxdt(mu0:mu0+17) = 2 ! Warning d2gxdt is here pure imaginary
2280 !$OMP END SINGLE
2281              if (parity) then
2282                do mua=1,6 ! strain
2283                  do mub=1,3 ! k
2284                    mu=mu0+(mub-1)+3*(mua-1)
2285 !$OMP SINGLE
2286                    buffer_i1 = zero
2287 !$OMP END SINGLE
2288 !$OMP DO &
2289 !$OMP REDUCTION(+:buffer_i1)
2290                    do ipw=1,npw
2291                      buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub)
2292                    end do
2293 !$OMP END DO
2294 !$OMP SINGLE
2295                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1
2296 !$OMP END SINGLE
2297                  end do
2298                end do
2299              else
2300                do mua=1,6 ! strain
2301                  do mub=1,3 ! k
2302                    mu=mu0+(mub-1)+3*(mua-1)
2303 !$OMP SINGLE
2304                    buffer_r1 = zero
2305 !$OMP END SINGLE
2306 !$OMP DO &
2307 !$OMP REDUCTION(+:buffer_r1)
2308                    do ipw=1,npw
2309                      buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub)
2310                    end do
2311 !$OMP END DO
2312 !$OMP SINGLE
2313                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1
2314 !$OMP END SINGLE
2315                  end do
2316                end do
2317              end if
2318            end if
2319          end do
2320        end if
2321 
2322 !      --------------------------------------------------------------------
2323 !      CHOICE= 6  --  SIGNS= 1
2324 !      Accumulate d2Gxdt --- 2nd derivative wrt atm. pos and strain --- for all dirs
2325 !        --------------------------------------------------------------------
2326        if ((signs==1).and.(choice_==6)) then
2327          mu0=1;if (choice_==6) mu0=37
2328          ABI_ALLOCATE(scalarr,(npw,4))
2329          ABI_ALLOCATE(scalari,(npw,4))
2330          do ilmn=1,nlmn
2331            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2332            scale=wt;if (il>1) scale=-scale
2333            scale2=pi*scale
2334            do mu=1,4
2335 !$OMP SINGLE
2336              do ipw=1,npw
2337                scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn)
2338                scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn)
2339              end do
2340 !$OMP END SINGLE
2341            end do
2342            if(cplex==2) then
2343              do mub=1,6
2344                nub1=alpha(mub);nub2=beta(mub)
2345                do mua=1,3
2346                  mu=mu0+(mua-1)+3*(mub-1)
2347                  buffer_r1 = zero ; buffer_i1 = zero
2348 !$OMP DO &
2349 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
2350                  do ipw=1,npw
2351                    buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + &
2352 &                   scalari(ipw,1+nub2)*kpg(ipw,nub1))
2353                    buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + &
2354 &                   scalarr(ipw,1+nub2)*kpg(ipw,nub1))
2355                  end do
2356 !$OMP END DO
2357                  if (parity) then
2358 !$OMP SINGLE
2359                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1
2360                    d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_r1
2361 !$OMP END SINGLE
2362                  else
2363 !$OMP SINGLE
2364                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2365                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1
2366 !$OMP END SINGLE
2367                  end if
2368                end do
2369              end do
2370            else
2371              if (parity) then
2372                do mub=1,6
2373                  nub1=alpha(mub);nub2=beta(mub)
2374                  do mua=1,3
2375                    mu=mu0+(mua-1)+3*(mub-1)
2376                    buffer_i1 = zero
2377 !$OMP DO &
2378 !$OMP REDUCTION(+:buffer_r1)
2379                    do ipw=1,npw
2380                      buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + &
2381 &                     scalari(ipw,1+nub2)*kpg(ipw,nub1))
2382                    end do
2383 !$OMP END DO
2384 !$OMP SINGLE
2385                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1
2386 !$OMP END SINGLE
2387                  end do
2388                end do
2389              else
2390                do mub=1,6
2391                  nub1=alpha(mub);nub2=beta(mub)
2392                  do mua=1,3
2393                    mu=mu0+(mua-1)+3*(mub-1)
2394                    buffer_r1 = zero
2395 !$OMP DO &
2396 !$OMP REDUCTION(+:buffer_r1)
2397                    do ipw=1,npw
2398                      buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + &
2399 &                     scalarr(ipw,1+nub2)*kpg(ipw,nub1))
2400                    end do
2401 !$OMP END DO
2402 !$OMP SINGLE
2403                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2404 !$OMP END SINGLE
2405                  end do
2406                end do
2407              end if
2408            end if
2409          end do
2410          ABI_DEALLOCATE(scalarr)
2411          ABI_DEALLOCATE(scalari)
2412        end if
2413 
2414 !      --------------------------------------------------------------------
2415 !      CHOICE= 6  --  SIGNS= 1
2416 !      Accumulate d2Gxdt --- 2nd derivative wrt 2 strains --- for all dirs
2417 !        --------------------------------------------------------------------
2418        if ((signs==1).and.(choice_==6)) then
2419          mu0=1
2420          ABI_ALLOCATE(scalarr,(npw,10))
2421          ABI_ALLOCATE(scalari,(npw,10))
2422          do ilmn=1,nlmn
2423            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2424            scale=wt;if (il>1) scale=-scale
2425            scale2=quarter*scale
2426            do mu=1,10
2427 !$OMP SINGLE
2428              do ipw=1,npw
2429                scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn)
2430                scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn)
2431              end do
2432 !$OMP END SINGLE
2433            end do
2434            if(cplex==2) then
2435              do mub=1,6
2436                nub1=alpha(mub);nub2=beta(mub)
2437                do mua=1,6
2438                  mu=mu0+(mua-1)+6*(mub-1)
2439                  nua1=alpha(mua);nua2=beta(mua)
2440                  gama=gamma(nub1,nua1)
2441                  gamb=gamma(nub2,nua1)
2442                  gamc=gamma(nub1,nua2)
2443                  gamd=gamma(nub2,nua2)
2444                  buffer_r1 = zero ; buffer_i1 = zero
2445 !$OMP DO &
2446 !$OMP REDUCTION(+:buffer_r1,buffer_i1), &
2447 !$OMP PRIVATE(kpga,kpgb,kpgc,kpgd)
2448                  do ipw=1,npw
2449                    kpga=kpg(ipw,nub1)
2450                    kpgb=kpg(ipw,nub2)
2451                    kpgc=kpg(ipw,nua1)
2452                    kpgd=kpg(ipw,nua2)
2453                    buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd &
2454 &                   + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc
2455                    buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd &
2456 &                   + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc
2457                  end do
2458 !$OMP END DO
2459                  if (parity) then
2460 !$OMP SINGLE
2461                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2462                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1
2463 !$OMP END SINGLE
2464                  else
2465 !$OMP SINGLE
2466                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
2467                    d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2468 !$OMP END SINGLE
2469                  end if
2470                end do
2471              end do
2472            else
2473              if (parity) then
2474                do mub=1,6
2475                  nub1=alpha(mub);nub2=beta(mub)
2476                  do mua=1,6
2477                    mu=mu0+(mua-1)+6*(mub-1)
2478                    nua1=alpha(mua);nua2=beta(mua)
2479                    gama=gamma(nub1,nua1)
2480                    gamb=gamma(nub2,nua1)
2481                    gamc=gamma(nub1,nua2)
2482                    gamd=gamma(nub2,nua2)
2483                    buffer_r1 = zero
2484 !$OMP DO &
2485 !$OMP REDUCTION(+:buffer_r1), &
2486 !$OMP PRIVATE(kpga,kpgb,kpgc,kpgd)
2487                    do ipw=1,npw
2488                      kpga=kpg(ipw,nub1)
2489                      kpgb=kpg(ipw,nub2)
2490                      kpgc=kpg(ipw,nua1)
2491                      kpgd=kpg(ipw,nua2)
2492                      buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd &
2493 &                     + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc
2494                    end do
2495 !$OMP END DO
2496 !$OMP SINGLE
2497                    d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1
2498 !$OMP END SINGLE
2499                  end do
2500                end do
2501              else
2502                do mub=1,6
2503                  nub1=alpha(mub);nub2=beta(mub)
2504                  do mua=1,6
2505                    mu=mu0+(mua-1)+6*(mub-1)
2506                    nua1=alpha(mua);nua2=beta(mua)
2507                    gama=gamma(nub1,nua1)
2508                    gamb=gamma(nub2,nua1)
2509                    gamc=gamma(nub1,nua2)
2510                    gamd=gamma(nub2,nua2)
2511                    buffer_i1 = zero
2512 !$OMP DO &
2513 !$OMP REDUCTION(+:buffer_r1)
2514                    do ipw=1,npw
2515                      kpga=kpg(ipw,nub1)
2516                      kpgb=kpg(ipw,nub2)
2517                      kpgc=kpg(ipw,nua1)
2518                      kpgd=kpg(ipw,nua2)
2519                      buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd &
2520 &                     + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc
2521                    end do
2522 !$OMP END DO
2523 !$OMP SINGLE
2524                    d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1
2525 !$OMP END SINGLE
2526                  end do
2527                end do
2528              end if
2529            end if
2530          end do
2531          ABI_DEALLOCATE(scalarr)
2532          ABI_DEALLOCATE(scalari)
2533        end if
2534 
2535 !      --------------------------------------------------------------------
2536 !      CHOICE= 8, 81  --  SIGNS= 1
2537 !      Accumulate d2Gxdt --- 2nd derivative wrt 2 k wave vect. --- for all dirs
2538 !      --------------------------------------------------------------------
2539        if ((signs==1).and.(choice_==8.or.choice_==81)) then
2540          mu0=1
2541          do ilmn=1,nlmn
2542            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2543            scale=wt;if (il>1) scale=-scale
2544            if (cplex==2) then
2545 !$OMP SINGLE
2546              buffer_ra = zero ; buffer_rb = zero
2547              buffer_rc = zero ; buffer_rd = zero
2548              buffer_re = zero ; buffer_rf = zero
2549              buffer_ia = zero ; buffer_ib = zero
2550              buffer_ic = zero ; buffer_id = zero
2551              buffer_ie = zero ; buffer_if = zero
2552 !$OMP END SINGLE
2553 !$OMP DO &
2554 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc), &
2555 !$OMP REDUCTION(+:buffer_rd,buffer_re,buffer_rf), &
2556 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic), &
2557 !$OMP REDUCTION(+:buffer_id,buffer_ie,buffer_if)
2558              do ipw=1,npw
2559                buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn)
2560                buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn)
2561                buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn)
2562                buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn)
2563                buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn)
2564                buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn)
2565                buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn)
2566                buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn)
2567                buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn)
2568                buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn)
2569                buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn)
2570                buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn)
2571              end do
2572 !$OMP END DO
2573 !$OMP SINGLE
2574              if (parity) then
2575                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_ra
2576                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_ia
2577                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb
2578                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_ib
2579                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc
2580                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_ic
2581                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd
2582                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_id
2583                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re
2584                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_ie
2585                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf
2586                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_if
2587              else
2588                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_ia
2589                d2gxdt(2,mu0  ,ilmn,ia,ispinor) = scale*buffer_ra
2590                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib
2591                d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb
2592                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic
2593                d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc
2594                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id
2595                d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd
2596                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie
2597                d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_re
2598                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if
2599                d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf
2600              end if
2601 !$OMP END SINGLE
2602            else
2603              if (parity) then
2604 !$OMP SINGLE
2605                buffer_ra = zero ; buffer_rb = zero
2606                buffer_rc = zero ; buffer_rd = zero
2607                buffer_re = zero ; buffer_rf = zero
2608 !$OMP END SINGLE
2609 !$OMP DO &
2610 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc), &
2611 !$OMP REDUCTION(+:buffer_rd,buffer_re,buffer_rf)
2612                do ipw=1,npw
2613                  buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn)
2614                  buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn)
2615                  buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn)
2616                  buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn)
2617                  buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn)
2618                  buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn)
2619                end do
2620 !$OMP END DO
2621 !$OMP SINGLE
2622                d2gxdt(1,mu0  ,ilmn,ia,ispinor) = scale*buffer_ra
2623                d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb
2624                d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc
2625                d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd
2626                d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re
2627                d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf
2628 !$OMP END SINGLE
2629              else
2630 !$OMP SINGLE
2631                buffer_ia = zero ; buffer_ib = zero
2632                buffer_ic = zero ; buffer_id = zero
2633                buffer_ie = zero ; buffer_if = zero
2634 !$OMP END SINGLE
2635 !$OMP DO &
2636 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic), &
2637 !$OMP REDUCTION(+:buffer_id,buffer_ie,buffer_if)
2638                do ipw=1,npw
2639                  buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn)
2640                  buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn)
2641                  buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn)
2642                  buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn)
2643                  buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn)
2644                  buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn)
2645                end do
2646 !$OMP END DO
2647 !$OMP SINGLE
2648                d2gxdt(1,mu0  ,ilmn,ia,ispinor) =-scale*buffer_ia
2649                d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib
2650                d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic
2651                d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id
2652                d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie
2653                d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if
2654 !$OMP END SINGLE
2655              end if
2656            end if
2657          end do
2658        end if
2659 
2660 !      --------------------------------------------------------------------
2661 !      CHOICE= 8, 81 --  SIGNS= 2
2662 !      Accumulate d2Gxdt --- 2nd derivative wrt k --- for direction IDIR (Voigt not.)
2663 !      --------------------------------------------------------------------
2664        if ((signs==2).and.(choice_==8.or.choice_==81)) then
2665          mu0=1
2666          !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
2667          mua=idir1(idir);mub=idir2(idir)
2668          !idir->(Voigt) (11->1,22->2,33->3,32->4,31->5,21->6)
2669          ffnl_dir(1)=gamma(mua,mub)
2670          do ilmn=1,nlmn
2671            il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
2672            scale=wt;if (il>1) scale=-scale
2673            if (cplex==2) then
2674 !$OMP SINGLE
2675              buffer_r1 = zero ; buffer_i1 = zero
2676 !$OMP END SINGLE
2677 !$OMP DO &
2678 !$OMP REDUCTION(+:buffer_r1,buffer_i1)
2679              do ipw=1,npw
2680                buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
2681                buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
2682              end do
2683 !$OMP END DO
2684 !$OMP SINGLE
2685              if (parity) then
2686                d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
2687                d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1
2688              else
2689                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
2690                d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1
2691              end if
2692 !$OMP END SINGLE
2693            else
2694              if (parity) then
2695 !$OMP SINGLE
2696                buffer_r1 = zero
2697 !$OMP END SINGLE
2698 !$OMP DO &
2699 !$OMP REDUCTION(+:buffer_r1)
2700                do ipw=1,npw
2701                  buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
2702                end do
2703 !$OMP END DO
2704 !$OMP SINGLE
2705                d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1
2706 !$OMP END SINGLE
2707              else
2708 !$OMP SINGLE
2709                buffer_i1 = zero
2710 !$OMP END SINGLE
2711 !$OMP DO &
2712 !$OMP REDUCTION(+:buffer_i1)
2713                do ipw=1,npw
2714                  buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn)
2715                end do
2716 !$OMP SINGLE
2717                d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1
2718 !$OMP END SINGLE
2719              end if
2720            end if
2721          end do
2722        end if
2723 
2724 !      --------------------------------------------------------------------
2725 !      END CHOICES
2726 !      --------------------------------------------------------------------
2727 
2728      end do ! End loop on atoms
2729 
2730 !    Deallocate temporary space
2731 !$OMP SECTIONS
2732 !$OMP SECTION
2733      ABI_DEALLOCATE(scali)
2734 !$OMP SECTION
2735      ABI_DEALLOCATE(scalr)
2736 !$OMP END SECTIONS
2737 
2738    end do !  End loop on spinorial components
2739 !$OMP END PARALLEL
2740 
2741 !  ==========================================================================
2742  end if
2743 
2744 !Has to reduce arrays in case of FFT parallelization
2745  if (mpi_enreg%nproc_fft>1) then
2746    call timab(48,1,tsec)
2747    if (choice>=0) then
2748      call xmpi_sum(gx,mpi_enreg%comm_fft,ierr)
2749    end if
2750    if (choice_>1) then
2751      if (ndgxdt>0) then
2752        call xmpi_sum(dgxdt,mpi_enreg%comm_fft,ierr)
2753      end if
2754      if (nd2gxdt>0) then
2755        call xmpi_sum(d2gxdt,mpi_enreg%comm_fft,ierr)
2756      end if
2757    end if
2758    call timab(48,2,tsec)
2759  end if
2760 
2761 end subroutine opernla_ylm