TABLE OF CONTENTS


ABINIT/m_opernl [ Modules ]

[ Top ] [ Modules ]

NAME

  m_opernl

FUNCTION

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, DRH)
  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

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_opernl
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_mkffkg, only : mkffkg, dfpt_mkffkg
34 
35  implicit none
36 
37  private

ABINIT/opernl2 [ Functions ]

[ Top ] [ Functions ]

NAME

 opernl2

FUNCTION

 Operate with the non-local part of the hamiltonian,
 either from reciprocal space to projected quantities (sign=1),
 or from projected quantities to reciprocal space (sign=-1)

INPUTS

  if(sign==-1 .and. (choice==2 .or choice==4 .or. choice==5))
   dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
    or with respect to ddk
  if(sign==-1 .and. choice==3)
   dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
  ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
  gmet(3,3)=metric tensor for G vecs (in bohr**-2)
  ia3=gives the number of the first atom in the subset presently treated
  idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
  ispinor=1 or 2, gives the spinorial component of ffnl to be used
  istwf_k=option parameter that describes the storage of wfs
  itypat = type of atom, needed for ffnl
  jproj(nlang)=number of projectors for each angular momentum
  kg_k(3,npw)=integer coords of planewaves in basis sphere
  kpg_k(npw,npkg)= (k+G) components and related data
  kpt(3)=real components of k point in terms of recip. translations
  lmnmax=max. number of (l,n) components over all type of psps
  matblk=dimension of the array ph3d
  mincat= maximum increment of atoms
  mlang1 = dimensions for dgxdis1
  mlang3 = one of the dimensions of the array gxa
  mlang4 = dimension for dgxds
  mlang5 = dimensions for dgxdis2
  mlang6 = dimension for d2gxds2
  mproj=maximum dimension for number of projection operators for each
    angular momentum for nonlocal pseudopotential
  ndgxdt=second dimension of dgxdt
  nffnl=second dimension of ffnl
  nincat = number of atoms in the subset here treated
  nkpg=second size of array kpg_k
  nlang = number of angular momenta to be treated = 1 + highest ang. mom.
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npw  = number of plane waves in reciprocal space
  ntypat = number of type of atoms, dimension needed for ffnl
  sign : if  1, go from reciprocal space to projected scalars,
         if -1, go from projected scalars to reciprocal space.
  if(sign==1), vect(2*npw)=starting vector in reciprocal space
  if(sign==-1) gxa(2,mlang3,nincat,mproj)=modified projected scalars;
   NOTE that metric contractions have already been performed on the
   arrays gxa if sign=-1
  ph3d(2,npw,matblk)=three-dimensional phase factors

OUTPUT

  if(sign==1)
   gxa(2,mlang3,mincat,mproj)= projected scalars
  if(sign==1 .and. (choice==2 .or choice==4 .or. choice==5 .or. choice==23))
   dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
    or with respect to ddk
  if(sign==1 .and. (choice==3 .or. choice==23))
   dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
  if(sign==1 .and. choice==6)
   dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars
    wrt coord. indexed for internal strain
   d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars
    wrt strain and coord
   d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars
    wrt strains
  if(sign==-1)
   vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.

NOTES

 Operate with the non-local part of the hamiltonian for one type of
 atom, and within this given type of atom, for a subset
 of at most nincat atoms.

 This routine basically replaces getgla (gxa here is the former gla),
 except for the calculation of <G|dVnl/dk|C> or strain gradients.

PARENTS

      nonlop_pl

CHILDREN

      mkffkg

SOURCE

136 subroutine opernl2(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
137 &  ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
138 &  jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
139 &  mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
140 &  ntypat,ph3d,sign,vect)
141 
142 
143 !This section has been created automatically by the script Abilint (TD).
144 !Do not modify the following lines by hand.
145 #undef ABI_FUNC
146 #define ABI_FUNC 'opernl2'
147 !End of the abilint section
148 
149  implicit none
150 
151 !Arguments ------------------------------------
152 !scalars
153  integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk
154  integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt
155  integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign
156 !arrays
157  integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
158  integer,intent(in) :: nloalg(3)
159  real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
160  real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk)
161  real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj)
162  real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
163  real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:)
164  real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj)
165  real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj)
166  real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj)
167 
168 !Local variables-------------------------------
169 !scalars
170  integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3
171  integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw
172  integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj,ntens
173  real(dp),parameter :: two_pi2=two_pi**2
174  character(len=500) :: message
175 !arrays
176  integer,allocatable :: parity(:)
177 !real(dp) :: tsec(2)
178  real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
179 
180 ! *************************************************************************
181 
182 !call wrtout(std_out,"in opernl2","COLL")
183 
184 !mblkpw sets the size of blocks of planewaves
185  mblkpw=NLO_MBLKPW
186 
187 !Get the actual maximum number of projectors
188  mmproj=maxval(indlmn(3,:,itypat))
189 
190 !Initialisation before blocking on the plane waves
191 
192  if (sign==1) then
193 !  Put projected scalars to zero
194    gxa(:,:,:,1:mmproj)=0.0d0
195    if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
196    if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
197    if (choice==6) then
198      dgxdis(:,:,:,1:mmproj)=0.0d0
199      d2gxdis(:,:,:,1:mmproj)=0.0d0
200      d2gxds2(:,:,:,1:mmproj)=0.0d0
201    end if
202  end if
203 
204 !Set up dimension of kpgx and allocate
205 !ntens sets the maximum number of independent tensor components
206 !over all allowed angular momenta; need 20 for spdf for tensors
207 !up to rank 3; to handle stress tensor, need up to rank 5
208  ntens=1
209  if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4
210  if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
211  if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
212  if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
213  if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
214  if(choice==6 .and. nlang>=3)ntens=84
215  if(choice==6 .and. nlang==4)ntens=120
216 
217 !Set up second dimension of ffkg array, and allocate
218  nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
219  do ilang=1,nlang
220 !  Get the number of projectors for that angular momentum
221    nproj=jproj(ilang)
222 !  If there is a non-local part, accumulate the number of vectors needed
223 !  The variables ilang below are the number of independent tensors of
224 !  various ranks, the variable names being more historical than logical.
225 !  ilang2=number of rank ilang-1
226 !  ilang3=number of rank ilang+1
227 !  ilang4=number of rank ilang
228 !  ilang5=number of rank ilang+2
229 !  ilang6=number of rank ilang+3
230    if(nproj>0)then
231      ilang2=(ilang*(ilang+1))/2
232      nffkge=nffkge+nproj*ilang2
233      if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
234      if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2
235      if(choice==3 .or. choice==6 .or. choice==23)then
236        ilang3=((ilang+2)*(ilang+3))/2
237        nffkgs=nffkgs+nproj*ilang3
238      end if
239      if(choice==6)then
240        ilang4=((ilang+1)*(ilang+2))/2
241        ilang5=((ilang+3)*(ilang+4))/2
242        ilang6=((ilang+4)*(ilang+5))/2
243        nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6)
244      end if
245    end if
246  end do
247  nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk
248 
249 !Loop on subsets of plane waves (blocking)
250 !Disabled by MG on Dec  6 2011, omp sections have to be tested, this coding causes a sigfault with nthreads==1
251 !Feb 16 2012: The code does not crash anymore but it's not efficient.
252 !
253 !!$OMP PARALLEL DEFAULT(PRIVATE) &
254 !!$OMP SHARED(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) &
255 !!$OMP SHARED(ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) &
256 !!$OMP SHARED(jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) &
257 !!$OMP SHARED(mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) &
258 !!$OMP SHARED(ntypat,ph3d,sign,vect) &
259 !!$OMP SHARED(mblkpw,nffkg,nffkgd,nffkge,nffkgs,ntens,mmproj)
260 
261  ABI_ALLOCATE(ffkg,(mblkpw,nffkg))
262  ABI_ALLOCATE(parity,(nffkg))
263  ABI_ALLOCATE(kpgx,(mblkpw,ntens))
264  ABI_ALLOCATE(scalars,(2,nffkg))
265  ABI_ALLOCATE(teffv,(2,mblkpw))
266 !!$OMP DO
267  do ipw1=1,npw,mblkpw
268 
269    ipw2=min(npw,ipw1+mblkpw-1)
270    nincpw=ipw2-ipw1+1
271 
272 !  call timab(74+choice,1,tsec)
273 
274 !  Initialize kpgx array related to tensors defined below
275    call mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,kg_k,&
276 &   kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
277 &   npw,ntens,ntypat,parity)
278 
279 !  call timab(74+choice,2,tsec)
280 
281 !  Now treat the different signs
282    if (sign==1) then
283 
284      do ia=1,nincat
285 
286 !      Compute the shift eventually needed to get the phases in ph3d
287        iaph3d=ia
288        if(nloalg(2)>0)iaph3d=ia+ia3-1
289 
290 !      ******* Entering the first time-consuming part of the routine *******
291 
292 !      Multiply by the phase factor
293 !      This allows to be left with only real operations,
294 !      that are performed in the most inner loops
295        ig=ipw1
296        do ipw=1,nincpw
297          teffv(1,ipw)=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
298          teffv(2,ipw)=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
299          ig=ig+1
300        end do
301 
302        do iffkg=1,nffkg
303          scalars(1,iffkg)=0.0d0
304          scalars(2,iffkg)=0.0d0
305          do ipw=1,nincpw
306            scalars(1,iffkg)=scalars(1,iffkg)+teffv(1,ipw)*ffkg(ipw,iffkg)
307            scalars(2,iffkg)=scalars(2,iffkg)+teffv(2,ipw)*ffkg(ipw,iffkg)
308          end do
309        end do
310 
311 !      ******* Leaving the critical part *********************************
312 
313 !      DEBUG
314 !      write(std_out,*)' opernl2 : scalars'
315 !      do iffkg=1,10
316 !      write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
317 !      end do
318 !      stop
319 !      ENDDEBUG
320 
321        if(istwf_k>=2)then
322 !        Impose parity of resulting scalar (this operation could be
323 !        replaced by direct saving of CPU time in the preceeding section)
324          do iffkg=1,nffkg
325            scalars(parity(iffkg),iffkg)=0.0d0
326          end do
327        end if
328 
329        iffkg=0
330        iffkgs=nffkge+nffkgd
331        iffkgs2=nffkge+nffkgs
332        iffkgk=nffkge*2
333        do ilang=1,nlang
334          nproj=jproj(ilang)
335          if(nproj>0)then
336 !          ilang2 is the number of independent tensor components
337 !          for symmetric tensor of rank ilang-1
338            ilang2=(ilang*(ilang+1))/2
339 
340 !          Loop over projectors
341            do iproj=1,nproj
342 !            Multiply by the k+G factors (tensors of various rank)
343              do ii=1,ilang2
344 !              Get the starting address for the relevant tensor
345                jj=ii+((ilang-1)*ilang*(ilang+1))/6
346                iffkg=iffkg+1
347 !              !$OMP CRITICAL (OPERNL2_1)
348                gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
349                gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg)
350 !              !$OMP END CRITICAL (OPERNL2_1)
351 !              Now, compute gradients, if needed.
352                if ((choice==2.or.choice==23) .and. ndgxdt==3) then
353                  do mu=1,3
354                    jffkg=nffkge+(iffkg-1)*3+mu
355 !                  Pay attention to the use of reals and imaginary parts here ...
356 !                  !$OMP CRITICAL (OPERNL2_2)
357                    dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
358                    dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
359 !                  !$OMP END CRITICAL (OPERNL2_2)
360                  end do
361                end if
362                if (choice==2 .and. ndgxdt==1) then
363                  jffkg=nffkge+iffkg
364 !                Pay attention to the use of reals and imaginary parts here ...
365 !                !$OMP CRITICAL (OPERNL2_3)
366                  dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg)
367                  dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg)
368 !                !$OMP END CRITICAL (OPERNL2_3)
369                end if
370                if (choice==4) then
371                  do mu=1,3
372                    jffkg=nffkge+(iffkg-1)*9+mu
373 !                  Pay attention to the use of reals and imaginary parts here ...
374 !                  !$OMP CRITICAL (OPERNL2_4)
375                    dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
376                    dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
377 !                  !$OMP END CRITICAL (OPERNL2_4)
378                  end do
379                  do mu=4,9
380                    jffkg=nffkge+(iffkg-1)*9+mu
381 !                  Pay attention to the use of reals and imaginary parts here ...
382 !                  Also, note the multiplication by (2 pi)**2
383 !                  !$OMP CRITICAL (OPERNL2_5)
384                    dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg)
385                    dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg)
386 !                  !$OMP END CRITICAL (OPERNL2_5)
387                  end do
388                end if
389 !              End loop on ii=1,ilang2
390              end do
391 
392              if (choice==3 .or. choice==6 .or. choice==23) then
393 !              Compute additional tensors related to strain gradients
394 !              ilang3 is number of unique tensor components of rank ilang+1
395                ilang3=((ilang+2)*(ilang+3))/2
396                jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
397 !              Compute strain gradient tensor components
398                do ii=1,ilang3
399 !                Note that iffkgs is also used by ddk and 2nd derivative parts
400                  iffkgs=iffkgs+1
401                  jj=ii+jjs
402 !                !$OMP CRITICAL (OPERNL2_6)
403                  dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
404                  dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs)
405 !                !$OMP END CRITICAL (OPERNL2_6)
406                end do
407              end if
408 
409              if (choice==6) then
410 !              Compute additional tensors related to strain 2nd derivatives
411 !              and internal strain derivatives
412 !              ilang6 is number of unique tensor components of rank ilang+3
413                ilang6=((ilang+4)*(ilang+5))/2
414                jjs=((ilang+3)*(ilang+4)*(ilang+5))/6
415 !              Compute strain gradient tensor components
416                do ii=1,ilang6
417                  iffkgs2=iffkgs2+1
418                  jj=ii+jjs
419 !                !$OMP CRITICAL (OPERNL2_6)
420                  d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2)
421                  d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2)
422 !                !$OMP END CRITICAL (OPERNL2_6)
423                end do
424 
425 !              ilang4 is number of unique tensor components of rank ilang
426                ilang4=((ilang+1)*(ilang+2))/2
427                jjs=((ilang)*(ilang+1)*(ilang+2))/6
428 !              Compute internal strain gradient tensor components
429                do ii=1,ilang4
430                  iffkgs2=iffkgs2+1
431                  jj=ii+jjs
432 !                !$OMP CRITICAL (OPERNL2_6)
433 !                Pay attention to the use of reals and imaginary parts here ...
434                  dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2)
435                  dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2)
436 !                !$OMP END CRITICAL (OPERNL2_6)
437                end do
438 
439 !              ilang5 is number of unique tensor components of rank ilang+2
440                ilang5=((ilang+3)*(ilang+4))/2
441                jjs=((ilang+2)*(ilang+3)*(ilang+4))/6
442 !              Compute internal strain gradient tensor components
443                do ii=1,ilang5
444                  iffkgs2=iffkgs2+1
445                  jj=ii+jjs
446 !                !$OMP CRITICAL (OPERNL2_6)
447 !                Pay attention to the use of reals and imaginary parts here ...
448                  d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2)
449                  d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2)
450 !                !$OMP END CRITICAL (OPERNL2_6)
451                end do
452              end if ! choice==6
453 
454              if (choice==5) then
455 !              Compute additional tensors related to ddk with ffnl(:,2,...)
456                ilangx=(ilang*(ilang+1))/2
457                jjs=((ilang-1)*ilang*(ilang+1))/6
458                do ii=1,ilangx
459 !                Note that iffkgs is also used by strain part
460                  iffkgs=iffkgs+1
461                  jj=ii+jjs
462 !                !$OMP CRITICAL (OPERNL2_7)
463                  dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
464                  dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs)
465 !                !$OMP END CRITICAL (OPERNL2_7)
466                end do
467 !              Compute additional tensors related to ddk with ffnl(:,1,...)
468                if(ilang>=2)then
469                  ilangx=((ilang-1)*ilang)/2
470                  jjs=((ilang-2)*(ilang-1)*ilang)/6
471                  do ii=1,ilangx
472                    iffkgk=iffkgk+1
473                    jj=ii+jjs
474 !                  !$OMP CRITICAL (OPERNL2_8)
475                    dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
476                    dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk)
477 !                  !$OMP END CRITICAL (OPERNL2_8)
478                  end do
479                end if
480              end if
481 
482 !            End projector loop
483            end do
484 
485 !          End condition of non-zero projectors
486          end if
487 
488 !        End angular momentum loop
489        end do
490 
491 !      End loop on atoms
492      end do
493 
494    else if (sign==-1) then
495 !    Application of non-local part from projected scalars
496 !    back to reciprocal space ...
497 !    [this section merely computes terms which add to <G|Vnl|C>;
498 !    nothing here is needed when various gradients are being computed]
499 
500 !    Loop on atoms
501      do ia=1,nincat
502 
503 !      Compute the shift eventually needed to get the phases in ph3d
504        iaph3d=ia
505        if(nloalg(2)>0)iaph3d=ia+ia3-1
506 
507 !      Transfer gxa (and eventually dgxdt) in scalars with different indexing
508        iffkg=0
509        iffkgk=nffkge*2
510        iffkgs=nffkge
511        do ilang=1,nlang
512          nproj=jproj(ilang)
513          if (nproj>0) then
514            ilang2=(ilang*(ilang+1))/2
515            ilang3=((ilang+2)*(ilang+3))/2
516            do iproj=1,nproj
517              do ii=1,ilang2
518                jj=ii+((ilang-1)*ilang*(ilang+1))/6
519                iffkg=iffkg+1
520                if(choice==1 .or. choice==3)then
521                  scalars(1,iffkg)=gxa(1,jj,ia,iproj)
522                  scalars(2,iffkg)=gxa(2,jj,ia,iproj)
523                else if (choice==2 .and. ndgxdt==1) then
524                  jffkg=nffkge+iffkg
525 !                Pay attention to the use of reals and imaginary parts here ...
526 !                Also, the gxa and dgxdt arrays are switched, in order
527 !                to give the correct combination when multiplying ffkg,
528 !                see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]]
529                  scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
530                  scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
531                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
532                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
533                else if (choice==5) then
534                  jffkg=nffkge+iffkg
535 !                The gxa and dgxdt arrays are switched, in order
536 !                to give the correct combination when multiplying ffkg,
537                  scalars(1,jffkg)= gxa(1,jj,ia,iproj)
538                  scalars(2,jffkg)= gxa(2,jj,ia,iproj)
539                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
540                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
541                end if
542              end do
543              if(choice==3) then
544                do ii=1,ilang3
545                  iffkgs=iffkgs+1
546                  jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
547                  scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
548                  scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
549                end do
550              end if
551              if(ilang>=2 .and. choice==5)then
552                do ii=1,((ilang-1)*ilang)/2
553                  jj=ii+((ilang-2)*(ilang-1)*ilang)/6
554                  iffkgk=iffkgk+1
555                  scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
556                  scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
557                end do
558              end if
559            end do
560          end if
561        end do
562 
563 !      DEBUG
564 !      if(choice==5)then
565 !      write(std_out,*)' opernl2 : write gxa(:,...) array '
566 !      do ii=1,10
567 !      write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
568 !      end do
569 !      write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
570 !      do ii=1,10
571 !      write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
572 !      end do
573 !      write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
574 !      do ii=1,4
575 !      write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
576 !      end do
577 !      end if
578 
579 !      do iffkg=1,nffkg
580 !      write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
581 !      end do
582 !      stop
583 !      ENDDEBUG
584 
585 !      ******* Entering the critical part *********************************
586 
587        do ipw=1,nincpw
588          teffv(1,ipw)=0.0d0
589          teffv(2,ipw)=0.0d0
590        end do
591        do iffkg=1,nffkg
592          do ipw=1,nincpw
593            teffv(1,ipw)=teffv(1,ipw)+ffkg(ipw,iffkg)*scalars(1,iffkg)
594            teffv(2,ipw)=teffv(2,ipw)+ffkg(ipw,iffkg)*scalars(2,iffkg)
595          end do
596        end do
597 !      Multiplication by the complex conjugate of the phase
598        ig=ipw1
599        do ipw=1,nincpw
600          vect(1,ig)=vect(1,ig)+teffv(1,ipw)*ph3d(1,ig,iaph3d)+teffv(2,ipw)*ph3d(2,ig,iaph3d)
601          vect(2,ig)=vect(2,ig)+teffv(2,ipw)*ph3d(1,ig,iaph3d)-teffv(1,ipw)*ph3d(2,ig,iaph3d)
602          ig=ig+1
603        end do
604 
605 !      ******* Leaving the critical part *********************************
606 
607 !      End loop on atoms
608      end do
609 
610 !    End sign==-1
611    else
612 
613 !    Problem: sign and/or choice do not make sense
614      write(message,'(a,2i10,a,a)')&
615 &     ' Input sign, choice=',sign,choice,ch10,&
616 &     ' Are not compatible or allowed. '
617      MSG_BUG(message)
618    end if
619 
620 !  End loop on blocks of planewaves
621  end do
622 !!$OMP END DO
623  ABI_DEALLOCATE(ffkg)
624  ABI_DEALLOCATE(kpgx)
625  ABI_DEALLOCATE(parity)
626  ABI_DEALLOCATE(scalars)
627  ABI_DEALLOCATE(teffv)
628 !!$OMP END PARALLEL
629 
630 
631 !DEBUG
632 !if(choice==5)then
633 !write(std_out,*)'opernl2 : write vect(2*npw)'
634 !do ipw=1,2
635 !write(std_out,*)ipw,vect(1:2,ipw)
636 !end do
637 !write(std_out,*)'opernl2 : write ph3d'
638 !do ipw=1,npw
639 !write(std_out,*)ipw,ph3d(1:2,ipw,1)
640 !end do
641 !write(std_out,*)' opernl2 : write gxa array '
642 !write(std_out,*)' ang mom ,ia '
643 !do iproj=1,mproj
644 !do ia=1,1
645 !do ii=1,3
646 !do ii=1,10
647 !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
648 !end do
649 !end do
650 !end do
651 !end if
652 !if(choice==5)then
653 !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
654 !do ii=1,10
655 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
656 !end do
657 !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
658 !do ii=1,4
659 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
660 !end do
661 !stop
662 !end if
663 !ENDDEBUG
664 
665 end subroutine opernl2

ABINIT/opernl3 [ Functions ]

[ Top ] [ Functions ]

NAME

 opernl3

FUNCTION

 Operate with the non-local part of the hamiltonian,
 either from reciprocal space to projected quantities (sign=1),
 or from projected quantities to reciprocal space (sign=-1)

INPUTS

  if(sign==-1 .and. (choice==2 .or choice==4 .or. choice==5))
   dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
    or with respect to ddk
  if(sign==-1 .and. choice==3)
   dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
  ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
  gmet(3,3)=metric tensor for G vecs (in bohr**-2)
  ia3=gives the number of the first atom in the subset presently treated
  idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
  ispinor=1 or 2, gives the spinorial component of ffnl to be used
  istwf_k=option parameter that describes the storage of wfs
  itypat = type of atom, needed for ffnl
  jproj(nlang)=number of projectors for each angular momentum
  kg_k(3,npw)=integer coords of planewaves in basis sphere
  kpg_k(npw,npkg)= (k+G) components and related data
  kpt(3)=real components of k point in terms of recip. translations
  lmnmax=max. number of (l,n) components over all type of psps
  matblk=dimension of the array ph3d
  mincat= maximum increment of atoms
  mlang1 = dimensions for dgxdis1
  mlang3 = one of the dimensions of the array gxa
  mlang4 = dimension for dgxds
  mlang5 = dimensions for dgxdis2
  mlang6 = dimension for d2gxds2
  mproj=maximum dimension for number of projection operators for each
    angular momentum for nonlocal pseudopotential
  ndgxdt=second dimension of dgxdt
  nffnl=third dimension of ffnl
  nincat = number of atoms in the subset here treated
  nkpg=second size of array kpg_k
  nlang = number of angular momenta to be treated = 1 + highest ang. mom.
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npw  = number of plane waves in reciprocal space
  ntypat = number of type of atoms, dimension needed for ffnl
  sign : if  1, go from reciprocal space to projected scalars,
         if -1, go from projected scalars to reciprocal space.
  if(sign==1),
   vect(2*npw)=starting vector in reciprocal space
  if(sign==-1)
   gxa(2,mlang3,nincat,mproj)=modified projected scalars;
   NOTE that metric contractions have already been performed on the
   arrays gxa if sign=-1
  ph3d(2,npw,matblk)=three-dimensional phase factors

OUTPUT

  if(sign==1)
   gxa(2,mlang3,mincat,mproj)= projected scalars
  if(sign==1 .and. (choice==2 .or choice==4 .or. choice==5 .or. choice==23))
   dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
    or with respect to ddk
  if(sign==1 .and. (choice==3 .or. choice==23))
   dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
  if(sign==1 .and. choice==6)
   dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars
    wrt coord. indexed for internal strain
   d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars
    wrt strain and coord
   d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars
    wrt strains
  if(sign==-1)
   vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.

NOTES

 Operate with the non-local part of the hamiltonian for one type of
 atom, and within this given type of atom, for a subset
 of at most nincat atoms.

 This routine basically replaces getgla (gxa here is the former gla),
 except for the calculation of <G|dVnl/dk|C> or strain gradients.

PARENTS

      nonlop_pl

CHILDREN

      dfpt_mkffkg

SOURCE

 757 subroutine opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
 758 &  ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
 759 &  jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 760 &  mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
 761 &  ntypat,ph3d,sign,vect)
 762 
 763 
 764 !This section has been created automatically by the script Abilint (TD).
 765 !Do not modify the following lines by hand.
 766 #undef ABI_FUNC
 767 #define ABI_FUNC 'opernl3'
 768 !End of the abilint section
 769 
 770  implicit none
 771 
 772 !Arguments ------------------------------------
 773 !scalars
 774  integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk
 775  integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt
 776  integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign
 777 !arrays
 778  integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
 779  integer,intent(in) :: nloalg(3)
 780  real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3)
 781  real(dp),intent(in) :: kpg_k(npw,nkpg),kpt(3),ph3d(2,npw,matblk)
 782  real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj)
 783  real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
 784  real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:)
 785  real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj)
 786  real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj)
 787  real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj)
 788 
 789 !Local variables-------------------------------
 790 !ntens sets the maximum number of independent tensor components
 791 !over all allowed angular momenta; need 20 for spdf for tensors
 792 !up to rank 3; to handle stress tensor, need up to rank 5
 793 !to handle strain 2DER, need up to rank 7
 794 !scalars
 795  integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3
 796  integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw
 797  integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj
 798  integer :: ntens
 799  real(dp) :: ai,ar
 800  real(dp),parameter :: two_pi2=two_pi**2
 801  character(len=500) :: message
 802 !arrays
 803  integer,allocatable :: parity(:)
 804 ! real(dp) :: tsec(2)
 805  real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
 806 
 807 ! *************************************************************************
 808 
 809 !call wrtout(std_out,"in opernl3","COLL")
 810 
 811 !mblkpw sets the size of blocks of planewaves
 812  mblkpw=NLO_MBLKPW
 813 
 814 !Get the actual maximum number of projectors
 815  mmproj=maxval(indlmn(3,:,itypat))
 816 
 817 !Initialisation before blocking on the plane waves
 818 
 819  if (sign==1) then
 820 !  Put projected scalars to zero
 821    gxa(:,:,:,1:mmproj)=0.0d0
 822    if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
 823    if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
 824    if (choice==6) then
 825      dgxdis(:,:,:,1:mmproj)=0.0d0
 826      d2gxdis(:,:,:,1:mmproj)=0.0d0
 827      d2gxds2(:,:,:,1:mmproj)=0.0d0
 828    end if
 829  end if
 830 
 831 !Set up dimension of kpgx and allocate
 832 !ntens sets the maximum number of independent tensor components
 833 !over all allowed angular momenta; need 20 for spdf for tensors
 834 !up to rank 3; to handle stress tensor, need up to rank 5
 835  ntens=1
 836  if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4
 837  if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
 838  if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
 839  if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
 840  if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
 841  if(choice==6 .and. nlang>=3)ntens=84
 842  if(choice==6 .and. nlang==4)ntens=120
 843 
 844 !Set up second dimension of ffkg array, and allocate
 845  nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
 846  do ilang=1,nlang
 847 !  Get the number of projectors for that angular momentum
 848    nproj=jproj(ilang)
 849 !  If there is a non-local part, accumulate the number of vectors needed
 850 !  The variables ilang below are the number of independent tensors of
 851 !  various ranks, the variable names being more historical than logical.
 852 !  ilang2=number of rank ilang-1
 853 !  ilang3=number of rank ilang+1
 854 !  ilang4=number of rank ilang
 855 !  ilang5=number of rank ilang+2
 856 !  ilang6=number of rank ilang+3
 857    if(nproj>0)then
 858      ilang2=(ilang*(ilang+1))/2
 859      nffkge=nffkge+nproj*ilang2
 860      if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
 861      if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2
 862      if(choice==3 .or. choice==6 .or. choice==23)then
 863        ilang3=((ilang+2)*(ilang+3))/2
 864        nffkgs=nffkgs+nproj*ilang3
 865      end if
 866      if(choice==6)then
 867        ilang4=((ilang+1)*(ilang+2))/2
 868        ilang5=((ilang+3)*(ilang+4))/2
 869        ilang6=((ilang+4)*(ilang+5))/2
 870        nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6)
 871      end if
 872    end if
 873  end do
 874  nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk
 875 
 876 !Disabled by MG on Dec  6 2011, omp sections have to be tested, this coding causes a
 877 !sigfault with nthreads==1
 878 !
 879 !Loop on subsets of plane waves (blocking)
 880 !!$OMP PARALLEL DEFAULT(PRIVATE) &
 881 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) &
 882 !!$OMP SHARED(istwf_k,itypat,jproj,kg_k,kpg_k,kpt,lmnmax,mblkpw,mproj) &
 883 !!$OMP SHARED(ndgxdt,nffkg,nffkgd,nffkge,nffkgs,nincat,nkpg,nlang) &
 884 !!$OMP SHARED(nloalg,ph3d,npw,ntens,ntypat,sign,vect)
 885 
 886  ABI_ALLOCATE(ffkg,(nffkg,mblkpw))
 887  ABI_ALLOCATE(parity,(nffkg))
 888  ABI_ALLOCATE(kpgx,(mblkpw,ntens))
 889  ABI_ALLOCATE(scalars,(2,nffkg))
 890  ABI_ALLOCATE(teffv,(2,mblkpw))
 891 !!$OMP DO
 892  do ipw1=1,npw,mblkpw
 893 
 894    ipw2=min(npw,ipw1+mblkpw-1)
 895    nincpw=ipw2-ipw1+1
 896 
 897 !  Initialize kpgx array related to tensors defined below
 898    call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
 899 &   kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
 900 &   npw,ntens,ntypat,parity)
 901 
 902 !  call timab(74+choice,2,tsec)
 903 
 904 !  Now treat the different signs
 905    if (sign==1) then
 906 
 907      do ia=1,nincat
 908 
 909 !      Compute the shift eventually needed to get the phases in ph3d
 910        iaph3d=ia
 911        if(nloalg(2)>0)iaph3d=ia+ia3-1
 912 
 913        do iffkg=1,nffkg
 914          scalars(1,iffkg)=0.0d0
 915          scalars(2,iffkg)=0.0d0
 916        end do
 917 
 918 !      ******* Entering the first time-consuming part of the routine *******
 919 
 920 !      Note : first multiply by the phase factor
 921 !      This allows to be left with only real operations afterwards.
 922        ig=ipw1
 923 
 924        do ipw=1,nincpw
 925          ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
 926          ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
 927          do iffkg=1,nffkg
 928            scalars(1,iffkg)=scalars(1,iffkg)+ar*ffkg(iffkg,ipw)
 929            scalars(2,iffkg)=scalars(2,iffkg)+ai*ffkg(iffkg,ipw)
 930          end do
 931          ig=ig+1
 932        end do
 933 
 934 !      ******* Leaving the critical part *********************************
 935 
 936 !      DEBUG
 937 !      write(std_out,*)' opernl2 : scalars'
 938 !      do iffkg=1,10
 939 !      write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
 940 !      end do
 941 !      stop
 942 !      ENDDEBUG
 943 
 944        if(istwf_k>=2)then
 945 !        Impose parity of resulting scalar (this operation could be
 946 !        replaced by direct saving of CPU time in the preceeding section)
 947          do iffkg=1,nffkg
 948            scalars(parity(iffkg),iffkg)=0.0d0
 949          end do
 950        end if
 951 
 952        iffkg=0
 953        iffkgs=nffkge+nffkgd
 954        iffkgs2=nffkge+nffkgs
 955        iffkgk=nffkge*2
 956        do ilang=1,nlang
 957          nproj=jproj(ilang)
 958          if(nproj>0)then
 959 !          ilang2 is the number of independent tensor components
 960 !          for symmetric tensor of rank ilang-1
 961            ilang2=(ilang*(ilang+1))/2
 962 
 963 !          Loop over projectors
 964            do iproj=1,nproj
 965 !            Multiply by the k+G factors (tensors of various rank)
 966              do ii=1,ilang2
 967 !              Get the starting address for the relevant tensor
 968                jj=ii+((ilang-1)*ilang*(ilang+1))/6
 969                iffkg=iffkg+1
 970 !              !$OMP CRITICAL (OPERNL3_1)
 971                gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
 972                gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg)
 973 !              !$OMP END CRITICAL (OPERNL3_1)
 974 !              Now, compute gradients, if needed.
 975                if ((choice==2.or.choice==23) .and. ndgxdt==3) then
 976                  do mu=1,3
 977                    jffkg=nffkge+(iffkg-1)*3+mu
 978 !                  Pay attention to the use of reals and imaginary parts here ...
 979 !                  !$OMP CRITICAL (OPERNL3_2)
 980                    dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
 981                    dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
 982 !                  !$OMP END CRITICAL (OPERNL3_2)
 983                  end do
 984                end if
 985                if (choice==2 .and. ndgxdt==1) then
 986                  jffkg=nffkge+iffkg
 987 !                Pay attention to the use of reals and imaginary parts here ...
 988 !                !$OMP CRITICAL (OPERNL3_3)
 989                  dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg)
 990                  dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg)
 991 !                !$OMP END CRITICAL (OPERNL3_3)
 992                end if
 993                if (choice==4) then
 994                  do mu=1,3
 995                    jffkg=nffkge+(iffkg-1)*9+mu
 996 !                  Pay attention to the use of reals and imaginary parts here ...
 997 !                  !$OMP CRITICAL (OPERNL3_4)
 998                    dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
 999                    dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
1000 !                  !$OMP END CRITICAL (OPERNL3_4)
1001                  end do
1002                  do mu=4,9
1003                    jffkg=nffkge+(iffkg-1)*9+mu
1004 !                  Pay attention to the use of reals and imaginary parts here ...
1005 !                  Also, note the multiplication by (2 pi)**2
1006 !                  !$OMP CRITICAL (OPERNL3_5)
1007                    dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg)
1008                    dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg)
1009 !                  !$OMP END CRITICAL (OPERNL3_5)
1010                  end do
1011                end if
1012 !              End loop on ii=1,ilang2
1013              end do
1014 
1015              if (choice==3 .or. choice==6 .or. choice==23) then
1016 !              Compute additional tensors related to strain gradients
1017 !              ilang3 is number of unique tensor components of rank ilang+1
1018                ilang3=((ilang+2)*(ilang+3))/2
1019                jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
1020 !              Compute strain gradient tensor components
1021                do ii=1,ilang3
1022 !                Note that iffkgs is also used by ddk and 2nd derivative parts
1023                  iffkgs=iffkgs+1
1024                  jj=ii+jjs
1025 !                !$OMP CRITICAL (OPERNL3_6)
1026                  dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
1027                  dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs)
1028 !                !$OMP END CRITICAL (OPERNL3_6)
1029                end do
1030              end if
1031 
1032              if (choice==6) then
1033 !              Compute additional tensors related to strain 2nd derivatives
1034 !              and internal strain derivatives
1035 !              ilang6 is number of unique tensor components of rank ilang+3
1036                ilang6=((ilang+4)*(ilang+5))/2
1037                jjs=((ilang+3)*(ilang+4)*(ilang+5))/6
1038 !              Compute strain gradient tensor components
1039                do ii=1,ilang6
1040 !                Note that iffkgs is also used by ddk part
1041                  iffkgs2=iffkgs2+1
1042                  jj=ii+jjs
1043 !                !$OMP CRITICAL (OPERNL3_6)
1044                  d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2)
1045                  d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2)
1046 !                !$OMP END CRITICAL (OPERNL3_6)
1047                end do
1048 
1049 !              ilang4 is number of unique tensor components of rank ilang
1050                ilang4=((ilang+1)*(ilang+2))/2
1051                jjs=((ilang)*(ilang+1)*(ilang+2))/6
1052 !              Compute internal strain gradient tensor components
1053                do ii=1,ilang4
1054                  iffkgs2=iffkgs2+1
1055                  jj=ii+jjs
1056 !                !$OMP CRITICAL (OPERNL3_6)
1057 !                Pay attention to the use of reals and imaginary parts here ...
1058                  dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2)
1059                  dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2)
1060 !                !$OMP END CRITICAL (OPERNL3_6)
1061                end do
1062 
1063 !              ilang5 is number of unique tensor components of rank ilang+2
1064                ilang5=((ilang+3)*(ilang+4))/2
1065                jjs=((ilang+2)*(ilang+3)*(ilang+4))/6
1066 !              Compute internal strain gradient tensor components
1067                do ii=1,ilang5
1068                  iffkgs2=iffkgs2+1
1069                  jj=ii+jjs
1070 !                !$OMP CRITICAL (OPERNL3_6)
1071 !                Pay attention to the use of reals and imaginary parts here ...
1072                  d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2)
1073                  d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2)
1074 !                !$OMP END CRITICAL (OPERNL3_6)
1075                end do
1076              end if ! choice==6
1077 
1078              if (choice==5) then
1079 !              Compute additional tensors related to ddk with ffnl(:,2,...)
1080                ilangx=(ilang*(ilang+1))/2
1081                jjs=((ilang-1)*ilang*(ilang+1))/6
1082                do ii=1,ilangx
1083 !                Note that iffkgs is also used by strain part
1084                  iffkgs=iffkgs+1
1085                  jj=ii+jjs
1086 !                !$OMP CRITICAL (OPERNL3_7)
1087                  dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
1088                  dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs)
1089 !                !$OMP END CRITICAL (OPERNL3_7)
1090                end do
1091 !              Compute additional tensors related to ddk with ffnl(:,1,...)
1092                if(ilang>=2)then
1093                  ilangx=((ilang-1)*ilang)/2
1094                  jjs=((ilang-2)*(ilang-1)*ilang)/6
1095                  do ii=1,ilangx
1096                    iffkgk=iffkgk+1
1097                    jj=ii+jjs
1098 !                  !$OMP CRITICAL (OPERNL3_8)
1099                    dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
1100                    dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk)
1101 !                  !$OMP END CRITICAL (OPERNL3_8)
1102                  end do
1103                end if
1104              end if
1105 
1106 !            End projector loop
1107            end do
1108 
1109 !          End condition of non-zero projectors
1110          end if
1111 
1112 !        End angular momentum loop
1113        end do
1114 
1115 !      End loop on atoms
1116      end do
1117 
1118    else if (sign==-1) then
1119 !    Application of non-local part from projected scalars
1120 !    back to reciprocal space ...
1121 !    [this section merely computes terms which add to <G|Vnl|C>;
1122 !    nothing here is needed when various gradients are being computed]
1123 
1124 !    Loop on atoms
1125      do ia=1,nincat
1126 
1127 !      Compute the shift eventually needed to get the phases in ph3d
1128        iaph3d=ia
1129        if(nloalg(2)>0)iaph3d=ia+ia3-1
1130 
1131 !      Transfer gxa (and eventually dgxdt) in scalars with different indexing
1132        iffkg=0
1133        iffkgk=nffkge*2
1134        iffkgs=nffkge
1135        do ilang=1,nlang
1136          nproj=jproj(ilang)
1137          if (nproj>0) then
1138            ilang2=(ilang*(ilang+1))/2
1139            ilang3=((ilang+2)*(ilang+3))/2
1140            do iproj=1,nproj
1141              do ii=1,ilang2
1142                jj=ii+((ilang-1)*ilang*(ilang+1))/6
1143                iffkg=iffkg+1
1144                if(choice==1 .or. choice==3)then
1145                  scalars(1,iffkg)=gxa(1,jj,ia,iproj)
1146                  scalars(2,iffkg)=gxa(2,jj,ia,iproj)
1147                else if (choice==2 .and. ndgxdt==1) then
1148                  jffkg=nffkge+iffkg
1149 !                Pay attention to the use of reals and imaginary parts here ...
1150 !                Also, the gxa and dgxdt arrays are switched, in order
1151 !                to give the correct combination when multiplying ffkg,
1152 !                see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]]
1153                  scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
1154                  scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
1155                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
1156                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
1157                else if (choice==5) then
1158                  jffkg=nffkge+iffkg
1159 !                The gxa and dgxdt arrays are switched, in order
1160 !                to give the correct combination when multiplying ffkg,
1161                  scalars(1,jffkg)= gxa(1,jj,ia,iproj)
1162                  scalars(2,jffkg)= gxa(2,jj,ia,iproj)
1163                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
1164                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
1165                end if
1166              end do
1167              if(choice==3) then
1168                do ii=1,ilang3
1169                  iffkgs=iffkgs+1
1170                  jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
1171                  scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
1172                  scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
1173                end do
1174              end if
1175              if(ilang>=2 .and. choice==5)then
1176                do ii=1,((ilang-1)*ilang)/2
1177                  jj=ii+((ilang-2)*(ilang-1)*ilang)/6
1178                  iffkgk=iffkgk+1
1179                  scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
1180                  scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
1181                end do
1182              end if
1183            end do
1184          end if
1185        end do
1186 
1187 !      DEBUG
1188 !      if(choice==5)then
1189 !      write(std_out,*)' opernl2 : write gxa(:,...) array '
1190 !      do ii=1,10
1191 !      write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
1192 !      end do
1193 !      write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
1194 !      do ii=1,10
1195 !      write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
1196 !      end do
1197 !      write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
1198 !      do ii=1,4
1199 !      write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
1200 !      end do
1201 !      end if
1202 
1203 !      do iffkg=1,nffkg
1204 !      write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
1205 !      end do
1206 !      stop
1207 !      ENDDEBUG
1208 
1209        ig=ipw1
1210 
1211 !      ******* Entering the critical part *********************************
1212 
1213        do ipw=1,nincpw
1214          ar=0.0d0
1215          ai=0.0d0
1216          do iffkg=1,nffkg
1217            ar=ar+ffkg(iffkg,ipw)*scalars(1,iffkg)
1218            ai=ai+ffkg(iffkg,ipw)*scalars(2,iffkg)
1219          end do
1220          vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
1221          vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
1222          ig=ig+1
1223        end do
1224 
1225 !      ******* Leaving the critical part *********************************
1226 
1227 !      End loop on atoms
1228      end do
1229 
1230 !    End sign==-1
1231    else
1232 
1233 !    Problem: sign and/or choice do not make sense
1234      write(message, '(a,2i10,a,a)' )&
1235 &     '  Input sign, choice=',sign,choice,ch10,&
1236 &     '  Are not compatible or allowed. '
1237      MSG_BUG(message)
1238    end if
1239 
1240 !  End loop on blocks of planewaves
1241  end do
1242 !!$OMP END DO
1243  ABI_DEALLOCATE(ffkg)
1244  ABI_DEALLOCATE(kpgx)
1245  ABI_DEALLOCATE(parity)
1246  ABI_DEALLOCATE(scalars)
1247  ABI_DEALLOCATE(teffv)
1248 !!$OMP END PARALLEL
1249 
1250 
1251 !DEBUG
1252 !if(choice==5)then
1253 !write(std_out,*)'opernl2 : write vect(2*npw)'
1254 !do ipw=1,2
1255 !write(std_out,*)ipw,vect(1:2,ipw)
1256 !end do
1257 !write(std_out,*)'opernl2 : write ph3d'
1258 !do ipw=1,npw
1259 !write(std_out,*)ipw,ph3d(1:2,ipw,1)
1260 !end do
1261 !write(std_out,*)' opernl2 : write gxa array '
1262 !write(std_out,*)' ang mom ,ia '
1263 !do iproj=1,mproj
1264 !do ia=1,1
1265 !do ii=1,3
1266 !do ii=1,10
1267 !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
1268 !end do
1269 !end do
1270 !end do
1271 !end if
1272 !if(choice==5)then
1273 !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
1274 !do ii=1,10
1275 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
1276 !end do
1277 !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
1278 !do ii=1,4
1279 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
1280 !end do
1281 !stop
1282 !end if
1283 !ENDDEBUG
1284 
1285 end subroutine opernl3

ABINIT/opernl4a [ Functions ]

[ Top ] [ Functions ]

NAME

 opernl4a

FUNCTION

 Operate with the non-local part of the hamiltonian,
 from reciprocal space to projected quantities
 (oprnl4b is from projected quantities to reciprocal space)

INPUTS

  ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
  ia3=gives the number of the first atom in the subset presently treated
  idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
  ispinor=1 or 2, gives the spinorial component of ffnl to be used
  istwf_k=option parameter that describes the storage of wfs
  itypat = type of atom, needed for ffnl
  jproj(nlang)=number of projectors for each angular momentum
  kg_k(3,npw)=integer coords of planewaves in basis sphere
  kpg_k(npw,npkg)= (k+G) components and related data
  kpt(3)=real components of k point in terms of recip. translations
  lmnmax=max. number of (l,n) components over all type of psps
  matblk=dimension of the array ph3d
  mincat= maximum increment of atoms
  mlang1 = dimensions for dgxdis1
  mlang3 = one of the dimensions of the array gxa
  mlang4 = dimension for dgxds
  mlang5 = dimensions for dgxdis2
  mlang6 = dimension for d2gxds2
  mproj=maximum dimension for number of projection operators for each
    angular momentum for nonlocal pseudopotential
  ndgxdt=second dimension of dgxdt
  nffnl=third dimension of ffnl
  nincat = number of atoms in the subset here treated
  nkpg=second size of array kpg_k
  nlang = number of angular momenta to be treated = 1 + highest ang. mom.
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npw  = number of plane waves in reciprocal space
  ntypat = number of type of atoms, dimension needed for ffnl
  vect(2*npw)=starting vector in reciprocal space
  ph3d(2,npw,matblk)=three-dimensional phase factors

OUTPUT

  gxa(2,mlang3,mincat,mproj)= projected scalars
  if(choice==2 .or choice==4 .or. choice==5 .or. choice==23)
   dgxdt(2,ndgxdt,mlang3,mincat,mproj)= gradients of projected scalars wrt coords
    or with respect to ddk
  if(choice==3 .or. choice==23)
   dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
  if(choice==6)
   dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars
    wrt coord. indexed for internal strain
   d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars
    wrt strain and coord
   d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars
    wrt strains

NOTES

 Operate with the non-local part of the hamiltonian for one type of
 atom, and within this given type of atom, for a subset
 of at most nincat atoms.
 This routine basically replaces getgla (gxa here is the former gla),
 except for the calculation of <G|dVnl/dk|C> or strain gradients.

 Present version decomposed according to iffkg

PARENTS

      nonlop_pl

CHILDREN

      dfpt_mkffkg

SOURCE

1363 subroutine opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
1364 &  ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
1365 &  jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1366 &  mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
1367 &  ntypat,ph3d,vect)
1368 
1369 
1370 !This section has been created automatically by the script Abilint (TD).
1371 !Do not modify the following lines by hand.
1372 #undef ABI_FUNC
1373 #define ABI_FUNC 'opernl4a'
1374 !End of the abilint section
1375 
1376  implicit none
1377 
1378 !Arguments ------------------------------------
1379 !scalars
1380  integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk
1381  integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt
1382  integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat
1383 !arrays
1384  integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
1385  integer,intent(in) :: nloalg(3)
1386  real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
1387  real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk),vect(:,:)
1388  real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj)
1389  real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj)
1390  real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj)
1391  real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) !vz_i
1392  real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) !vz_i
1393  real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj) !vz_i
1394 
1395 !Local variables-------------------------------
1396 !scalars
1397  integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2
1398  integer :: ilang3,ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs
1399  integer :: jump,mblkpw,mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2
1400  integer :: nincpw,nproj,ntens,start
1401  real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8
1402  real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8
1403  real(dp),parameter :: two_pi2=two_pi*two_pi
1404 !arrays
1405  integer,allocatable :: parity(:)
1406  real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
1407 
1408 ! *************************************************************************
1409 
1410 !call wrtout(std_out,"in opernl4a","COLL")
1411 
1412 !mblkpw sets the size of blocks of planewaves
1413  mblkpw=NLO_MBLKPW
1414 
1415 !jump governs, in fine, the use of registers in the most cpu
1416 !time consuming part of the routine. Until now, jump=8 is the maximal value.
1417 !The optimal value will be machine-dependent !
1418  jump=4
1419 
1420 !Get the actual maximum number of projectors
1421  mmproj=maxval(indlmn(3,:,itypat))
1422 
1423 !Initialisation before blocking on the plane waves
1424 
1425 !Put projected scalars to zero
1426  gxa(:,:,:,1:mmproj)=0.0d0
1427  if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
1428  if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
1429  if (choice==6) then
1430    dgxdis(:,:,:,1:mmproj)=0.0d0
1431    d2gxdis(:,:,:,1:mmproj)=0.0d0
1432    d2gxds2(:,:,:,1:mmproj)=0.0d0
1433  end if
1434 
1435 !Set up dimension of kpgx and allocate
1436 !ntens sets the maximum number of independent tensor components
1437 !over all allowed angular momenta; need 20 for spdf for tensors
1438 !up to rank 3; to handle stress tensor, need up to rank 5
1439  ntens=1
1440  if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23)ntens=4
1441  if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
1442  if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
1443  if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
1444  if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
1445  if(choice==6 .and. nlang>=3)ntens=84
1446  if(choice==6 .and. nlang==4)ntens=120
1447 
1448 !Set up second dimension of ffkg array, and allocate
1449  nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
1450  do ilang=1,nlang
1451 !  Get the number of projectors for that angular momentum
1452    nproj=jproj(ilang)
1453 !  If there is a non-local part, accumulate the number of vectors needed
1454 !  The variables ilang below are the number of independent tensors of
1455 !  various ranks, the variable names being more historical than logical.
1456 !  ilang2=number of rank ilang-1
1457 !  ilang3=number of rank ilang+1
1458 !  ilang4=number of rank ilang
1459 !  ilang5=number of rank ilang+2
1460 !  ilang6=number of rank ilang+3
1461    if(nproj>0)then
1462      ilang2=(ilang*(ilang+1))/2
1463      nffkge=nffkge+nproj*ilang2
1464      if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
1465      if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2
1466      if(choice==3 .or. choice==6 .or. choice==23)then
1467        ilang3=((ilang+2)*(ilang+3))/2
1468        nffkgs=nffkgs+nproj*ilang3
1469      end if
1470      if(choice==6)then
1471        ilang4=((ilang+1)*(ilang+2))/2
1472        ilang5=((ilang+3)*(ilang+4))/2
1473        ilang6=((ilang+4)*(ilang+5))/2
1474        nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6)
1475      end if
1476    end if
1477  end do
1478  nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk
1479 
1480 !DEBUG
1481 !write(std_out,*)' jproj(1:nlang)',jproj(1:nlang)
1482 !write(std_out,*)' nffkg,nffkge,nffkgd,nffkgs,nffkgk',nffkg,nffkge,nffkgd,nffkgs,nffkgk
1483 !ENDDEBUG
1484 
1485 !Loop on subsets of plane waves (blocking)
1486 
1487 !Disabled by MG on Dec  6 2011, omp sections have to be tested, this coding causes a
1488 !sigfault with nthreads==1
1489 !Feb 16 2012: The code does not crash anymore but it's not efficient.
1490 !
1491 !!$OMP PARALLEL DEFAULT(PRIVATE) &
1492 !!$OMP SHARED (choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) &
1493 !!$OMP SHARED (ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) &
1494 !!$OMP SHARED (jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) &
1495 !!$OMP SHARED (mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) &
1496 !!$OMP SHARED (ntypat,ph3d,vect) &
1497 !!$OMP SHARED (mblkpw,jump,nffkgd,nffkg,nffkge,nffkgs,ntens)
1498 
1499  ABI_ALLOCATE(ffkg,(nffkg,mblkpw))
1500  ABI_ALLOCATE(parity,(nffkg))
1501  ABI_ALLOCATE(kpgx,(mblkpw,ntens))
1502  ABI_ALLOCATE(scalars,(2,nffkg))
1503  ABI_ALLOCATE(teffv,(2,mblkpw))
1504 
1505 !!$OMP DO
1506  do ipw1=1,npw,mblkpw
1507 
1508    ipw2=min(npw,ipw1+mblkpw-1)
1509    nincpw=ipw2-ipw1+1
1510 
1511 !  Initialize kpgx array related to tensors defined below
1512    call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
1513 &   kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
1514 &   npw,ntens,ntypat,parity)
1515 
1516    do ia=1,nincat
1517 
1518 !    Compute the shift eventually needed to get the phases in ph3d
1519      iaph3d=ia
1520      if(nloalg(2)>0)iaph3d=ia+ia3-1
1521 
1522      do iffkg=1,nffkg
1523        scalars(1,iffkg)=0.0d0 ; scalars(2,iffkg)=0.0d0
1524      end do
1525 
1526 !    DEBUG
1527 !    write(std_out,*)'opernl4, before first time-consuming'
1528 !    write(std_out,*)'opernl4 : nffkg,nincpw=',nffkg,nincpw
1529 !    write(std_out,*)'ig,ipw,ffkg(1:4),vect(1:2)'
1530 !    ig=ipw1
1531 !    do ipw=1,nincpw
1532 !    write(std_out,'(2i4,13es11.3)' )ig,ipw,ffkg(1:min(9,nffkg),ipw),vect(1:2,ipw),ph3d(1:2,ipw,iaph3d)
1533 !    ig=ig+1
1534 !    end do
1535 !    stop
1536 !    ENDDEBUG
1537 
1538 !    ******* Entering the first time-consuming part of the routine *******
1539 
1540 
1541 !    First, treat small nffkg; send treat the initial phase of big
1542 !    nffkg; finally treat the loop needed for big nffkg
1543 
1544 !    In the loops, first multiply by the phase factor.
1545 !    This allows to be left with only real operations afterwards.
1546 
1547 !    For the time being, the maximal jump allowed is 8.
1548 
1549 !    1) Here, treat small nffkg
1550      if(nffkg<=jump)then
1551 
1552        select case(nffkg)
1553 
1554        case(1)
1555 
1556          scr1=0.0d0 ; sci1=0.0d0
1557          ig=ipw1
1558          do ipw=1,nincpw
1559            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1560            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1561            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1562            ig=ig+1
1563          end do
1564          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1565 
1566        case(2)
1567 
1568          ig=ipw1
1569          scr1=0.0d0 ; sci1=0.0d0
1570          scr2=0.0d0 ; sci2=0.0d0
1571          do ipw=1,nincpw
1572            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1573            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1574            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1575            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1576            ig=ig+1
1577          end do
1578          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1579          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1580 
1581        case(3)
1582 
1583          ig=ipw1
1584          scr1=0.0d0 ; sci1=0.0d0
1585          scr2=0.0d0 ; sci2=0.0d0
1586          scr3=0.0d0 ; sci3=0.0d0
1587          do ipw=1,nincpw
1588            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1589            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1590            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1591            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1592            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1593            ig=ig+1
1594          end do
1595          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1596          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1597          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1598 
1599        case(4)
1600 
1601          ig=ipw1
1602          scr1=0.0d0 ; sci1=0.0d0
1603          scr2=0.0d0 ; sci2=0.0d0
1604          scr3=0.0d0 ; sci3=0.0d0
1605          scr4=0.0d0 ; sci4=0.0d0
1606          do ipw=1,nincpw
1607            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1608            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1609            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1610            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1611            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1612            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1613            ig=ig+1
1614          end do
1615          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1616          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1617          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1618          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1619 
1620        case(5)
1621 
1622          ig=ipw1
1623          scr1=0.0d0 ; sci1=0.0d0
1624          scr2=0.0d0 ; sci2=0.0d0
1625          scr3=0.0d0 ; sci3=0.0d0
1626          scr4=0.0d0 ; sci4=0.0d0
1627          scr5=0.0d0 ; sci5=0.0d0
1628          do ipw=1,nincpw
1629            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1630            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1631            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1632            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1633            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1634            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1635            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1636            ig=ig+1
1637          end do
1638          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1639          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1640          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1641          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1642          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1643 
1644        case(6)
1645 
1646          ig=ipw1
1647          scr1=0.0d0 ; sci1=0.0d0
1648          scr2=0.0d0 ; sci2=0.0d0
1649          scr3=0.0d0 ; sci3=0.0d0
1650          scr4=0.0d0 ; sci4=0.0d0
1651          scr5=0.0d0 ; sci5=0.0d0
1652          scr6=0.0d0 ; sci6=0.0d0
1653          do ipw=1,nincpw
1654            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1655            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1656            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1657            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1658            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1659            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1660            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1661            scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1662            ig=ig+1
1663          end do
1664          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1665          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1666          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1667          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1668          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1669          scalars(1,6)=scr6 ; scalars(2,6)=sci6
1670 
1671        case(7)
1672 
1673          ig=ipw1
1674          scr1=0.0d0 ; sci1=0.0d0
1675          scr2=0.0d0 ; sci2=0.0d0
1676          scr3=0.0d0 ; sci3=0.0d0
1677          scr4=0.0d0 ; sci4=0.0d0
1678          scr5=0.0d0 ; sci5=0.0d0
1679          scr6=0.0d0 ; sci6=0.0d0
1680          scr7=0.0d0 ; sci7=0.0d0
1681          do ipw=1,nincpw
1682            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1683            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1684            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1685            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1686            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1687            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1688            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1689            scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1690            scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1691            ig=ig+1
1692          end do
1693          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1694          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1695          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1696          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1697          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1698          scalars(1,6)=scr6 ; scalars(2,6)=sci6
1699          scalars(1,7)=scr7 ; scalars(2,7)=sci7
1700 
1701        case(8)
1702 
1703          ig=ipw1
1704          scr1=0.0d0 ; sci1=0.0d0
1705          scr2=0.0d0 ; sci2=0.0d0
1706          scr3=0.0d0 ; sci3=0.0d0
1707          scr4=0.0d0 ; sci4=0.0d0
1708          scr5=0.0d0 ; sci5=0.0d0
1709          scr6=0.0d0 ; sci6=0.0d0
1710          scr7=0.0d0 ; sci7=0.0d0
1711          scr8=0.0d0 ; sci8=0.0d0
1712          do ipw=1,nincpw
1713            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1714            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1715            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1716            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1717            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1718            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1719            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1720            scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1721            scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1722            scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw)
1723            ig=ig+1
1724          end do
1725          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1726          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1727          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1728          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1729          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1730          scalars(1,6)=scr6 ; scalars(2,6)=sci6
1731          scalars(1,7)=scr7 ; scalars(2,7)=sci7
1732          scalars(1,8)=scr8 ; scalars(2,8)=sci8
1733 
1734        end select
1735 
1736      else
1737 !      Now treat big nffkg
1738 
1739 !      2) Here, initialize big nffkg. The only difference with the
1740 !      preceeding case is that the intermediate results are stored.
1741 
1742        select case(jump)
1743 
1744        case(1)
1745 
1746          scr1=0.0d0 ; sci1=0.0d0
1747          ig=ipw1
1748          do ipw=1,nincpw
1749            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1750            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1751            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1752            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1753            ig=ig+1
1754          end do
1755          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1756 
1757        case(2)
1758 
1759          ig=ipw1
1760          scr1=0.0d0 ; sci1=0.0d0
1761          scr2=0.0d0 ; sci2=0.0d0
1762          do ipw=1,nincpw
1763            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1764            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1765            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1766            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1767            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1768            ig=ig+1
1769          end do
1770          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1771          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1772 
1773        case(3)
1774 
1775          ig=ipw1
1776          scr1=0.0d0 ; sci1=0.0d0
1777          scr2=0.0d0 ; sci2=0.0d0
1778          scr3=0.0d0 ; sci3=0.0d0
1779          do ipw=1,nincpw
1780            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1781            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1782            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1783            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1784            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1785            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1786            ig=ig+1
1787          end do
1788          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1789          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1790          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1791 
1792        case(4)
1793 
1794          ig=ipw1
1795          scr1=0.0d0 ; sci1=0.0d0
1796          scr2=0.0d0 ; sci2=0.0d0
1797          scr3=0.0d0 ; sci3=0.0d0
1798          scr4=0.0d0 ; sci4=0.0d0
1799          do ipw=1,nincpw
1800            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1801            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1802            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1803            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1804            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1805            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1806            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1807            ig=ig+1
1808          end do
1809          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1810          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1811          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1812          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1813 
1814        case(5)
1815 
1816          ig=ipw1
1817          scr1=0.0d0 ; sci1=0.0d0
1818          scr2=0.0d0 ; sci2=0.0d0
1819          scr3=0.0d0 ; sci3=0.0d0
1820          scr4=0.0d0 ; sci4=0.0d0
1821          scr5=0.0d0 ; sci5=0.0d0
1822          do ipw=1,nincpw
1823            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1824            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1825            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1826            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1827            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1828            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1829            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1830            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1831            ig=ig+1
1832          end do
1833          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1834          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1835          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1836          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1837          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1838 
1839        case(6)
1840 
1841          ig=ipw1
1842          scr1=0.0d0 ; sci1=0.0d0
1843          scr2=0.0d0 ; sci2=0.0d0
1844          scr3=0.0d0 ; sci3=0.0d0
1845          scr4=0.0d0 ; sci4=0.0d0
1846          scr5=0.0d0 ; sci5=0.0d0
1847          scr6=0.0d0 ; sci6=0.0d0
1848          do ipw=1,nincpw
1849            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1850            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1851            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1852            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1853            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1854            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1855            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1856            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1857            scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1858            ig=ig+1
1859          end do
1860          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1861          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1862          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1863          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1864          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1865          scalars(1,6)=scr6 ; scalars(2,6)=sci6
1866 
1867        case(7)
1868 
1869          ig=ipw1
1870          scr1=0.0d0 ; sci1=0.0d0
1871          scr2=0.0d0 ; sci2=0.0d0
1872          scr3=0.0d0 ; sci3=0.0d0
1873          scr4=0.0d0 ; sci4=0.0d0
1874          scr5=0.0d0 ; sci5=0.0d0
1875          scr6=0.0d0 ; sci6=0.0d0
1876          scr7=0.0d0 ; sci7=0.0d0
1877          do ipw=1,nincpw
1878            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1879            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1880            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1881            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1882            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1883            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1884            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1885            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1886            scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1887            scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1888            ig=ig+1
1889          end do
1890          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1891          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1892          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1893          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1894          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1895          scalars(1,6)=scr6 ; scalars(2,6)=sci6
1896          scalars(1,7)=scr7 ; scalars(2,7)=sci7
1897 
1898        case(8)
1899 
1900          ig=ipw1
1901          scr1=0.0d0 ; sci1=0.0d0
1902          scr2=0.0d0 ; sci2=0.0d0
1903          scr3=0.0d0 ; sci3=0.0d0
1904          scr4=0.0d0 ; sci4=0.0d0
1905          scr5=0.0d0 ; sci5=0.0d0
1906          scr6=0.0d0 ; sci6=0.0d0
1907          scr7=0.0d0 ; sci7=0.0d0
1908          scr8=0.0d0 ; sci8=0.0d0
1909          do ipw=1,nincpw
1910            ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1911            ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1912            teffv(1,ipw)=ar          ; teffv(2,ipw)=ai
1913            scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1914            scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1915            scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1916            scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1917            scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1918            scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1919            scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1920            scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw)
1921            ig=ig+1
1922          end do
1923          scalars(1,1)=scr1 ; scalars(2,1)=sci1
1924          scalars(1,2)=scr2 ; scalars(2,2)=sci2
1925          scalars(1,3)=scr3 ; scalars(2,3)=sci3
1926          scalars(1,4)=scr4 ; scalars(2,4)=sci4
1927          scalars(1,5)=scr5 ; scalars(2,5)=sci5
1928          scalars(1,6)=scr6 ; scalars(2,6)=sci6
1929          scalars(1,7)=scr7 ; scalars(2,7)=sci7
1930          scalars(1,8)=scr8 ; scalars(2,8)=sci8
1931 
1932        end select
1933 
1934 !      3) Here, do-loop for big nffkg.
1935 
1936        do start=1+jump,nffkg,jump
1937          chunk=min(jump,nffkg-start+1)
1938 
1939          select case(chunk)
1940 
1941          case(1)
1942 
1943            scr1=0.0d0 ; sci1=0.0d0
1944            do ipw=1,nincpw
1945              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
1946              scr1=scr1+ar*ffkg(start,ipw) ; sci1=sci1+ai*ffkg(start,ipw)
1947            end do
1948            scalars(1,start)=scr1 ; scalars(2,start)=sci1
1949 
1950          case(2)
1951 
1952            scr1=0.0d0 ; sci1=0.0d0
1953            scr2=0.0d0 ; sci2=0.0d0
1954            do ipw=1,nincpw
1955              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
1956              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
1957              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1958            end do
1959            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
1960            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1961 
1962          case(3)
1963 
1964            scr1=0.0d0 ; sci1=0.0d0
1965            scr2=0.0d0 ; sci2=0.0d0
1966            scr3=0.0d0 ; sci3=0.0d0
1967            do ipw=1,nincpw
1968              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
1969              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
1970              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1971              scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1972            end do
1973            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
1974            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1975            scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1976 
1977          case(4)
1978 
1979            scr1=0.0d0 ; sci1=0.0d0
1980            scr2=0.0d0 ; sci2=0.0d0
1981            scr3=0.0d0 ; sci3=0.0d0
1982            scr4=0.0d0 ; sci4=0.0d0
1983            do ipw=1,nincpw
1984              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
1985              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
1986              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1987              scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1988              scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
1989            end do
1990            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
1991            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1992            scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1993            scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
1994 
1995          case(5)
1996 
1997            scr1=0.0d0 ; sci1=0.0d0
1998            scr2=0.0d0 ; sci2=0.0d0
1999            scr3=0.0d0 ; sci3=0.0d0
2000            scr4=0.0d0 ; sci4=0.0d0
2001            scr5=0.0d0 ; sci5=0.0d0
2002            do ipw=1,nincpw
2003              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
2004              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
2005              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
2006              scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
2007              scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
2008              scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
2009            end do
2010            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
2011            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
2012            scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
2013            scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
2014            scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
2015 
2016          case(6)
2017 
2018            scr1=0.0d0 ; sci1=0.0d0
2019            scr2=0.0d0 ; sci2=0.0d0
2020            scr3=0.0d0 ; sci3=0.0d0
2021            scr4=0.0d0 ; sci4=0.0d0
2022            scr5=0.0d0 ; sci5=0.0d0
2023            scr6=0.0d0 ; sci6=0.0d0
2024            do ipw=1,nincpw
2025              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
2026              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
2027              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
2028              scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
2029              scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
2030              scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
2031              scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw)
2032            end do
2033            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
2034            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
2035            scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
2036            scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
2037            scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
2038            scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6
2039 
2040          case(7)
2041 
2042            scr1=0.0d0 ; sci1=0.0d0
2043            scr2=0.0d0 ; sci2=0.0d0
2044            scr3=0.0d0 ; sci3=0.0d0
2045            scr4=0.0d0 ; sci4=0.0d0
2046            scr5=0.0d0 ; sci5=0.0d0
2047            scr6=0.0d0 ; sci6=0.0d0
2048            scr7=0.0d0 ; sci7=0.0d0
2049            do ipw=1,nincpw
2050              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
2051              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
2052              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
2053              scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
2054              scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
2055              scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
2056              scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw)
2057              scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw)
2058            end do
2059            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
2060            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
2061            scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
2062            scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
2063            scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
2064            scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6
2065            scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7
2066 
2067          case(8)
2068 
2069            scr1=0.0d0 ; sci1=0.0d0
2070            scr2=0.0d0 ; sci2=0.0d0
2071            scr3=0.0d0 ; sci3=0.0d0
2072            scr4=0.0d0 ; sci4=0.0d0
2073            scr5=0.0d0 ; sci5=0.0d0
2074            scr6=0.0d0 ; sci6=0.0d0
2075            scr7=0.0d0 ; sci7=0.0d0
2076            scr8=0.0d0 ; sci8=0.0d0
2077            do ipw=1,nincpw
2078              ar=teffv(1,ipw)                ; ai=teffv(2,ipw)
2079              scr1=scr1+ar*ffkg(start  ,ipw) ; sci1=sci1+ai*ffkg(start  ,ipw)
2080              scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
2081              scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
2082              scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
2083              scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
2084              scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw)
2085              scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw)
2086              scr8=scr8+ar*ffkg(start+7,ipw) ; sci8=sci8+ai*ffkg(start+7,ipw)
2087            end do
2088            scalars(1,start  )=scr1 ; scalars(2,start  )=sci1
2089            scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
2090            scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
2091            scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
2092            scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
2093            scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6
2094            scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7
2095            scalars(1,start+7)=scr8 ; scalars(2,start+7)=sci8
2096 
2097          end select
2098 
2099 !        End loop on start
2100        end do
2101 
2102 !      End if statement for small or big nffkg
2103      end if
2104 
2105 !    ******* Leaving the critical part *********************************
2106 
2107 !    DEBUG
2108 !    write(std_out,*)' opernl4a, write scalars '
2109 !    do iffkg=1,nffkg
2110 !    write(std_out,*)iffkg,scalars(1:2,iffkg)
2111 !    end do
2112 !    ENDDEBUG
2113 
2114      if(istwf_k>=2)then
2115 !      Impose parity of resulting scalar (this operation could be
2116 !      replaced by direct saving of CPU time in the preceeding section)
2117        do iffkg=1,nffkg
2118          scalars(parity(iffkg),iffkg)=0.0d0
2119        end do
2120      end if
2121 
2122      iffkg=0 ; iffkgs=nffkge+nffkgd ; iffkgk=nffkge*2
2123      iffkgs2=nffkge+nffkgs
2124      do ilang=1,nlang
2125        nproj=jproj(ilang)
2126        if(nproj>0)then
2127 !        ilang2 is the number of independent tensor components
2128 !        for symmetric tensor of rank ilang-1
2129          ilang2=(ilang*(ilang+1))/2
2130 
2131 !        Loop over projectors
2132          do iproj=1,nproj
2133 !          Multiply by the k+G factors (tensors of various rank)
2134            do ii=1,ilang2
2135 !            Get the starting address for the relevant tensor
2136              jj=ii+((ilang-1)*ilang*(ilang+1))/6
2137              iffkg=iffkg+1
2138 !            !$OMP CRITICAL (OPERNL4a_1)
2139              gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
2140              gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg)
2141 !            !$OMP END CRITICAL (OPERNL4a_1)
2142 !            Now, compute gradients, if needed.
2143              if ((choice==2.or.choice==23) .and. ndgxdt==3) then
2144                do mu=1,3
2145                  jffkg=nffkge+(iffkg-1)*3+mu
2146 !                Pay attention to the use of reals and imaginary parts here ...
2147 !                !$OMP CRITICAL (OPERNL4a_2)
2148                  dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2149                  dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
2150 !                !$OMP END CRITICAL (OPERNL4a_2)
2151                end do
2152              end if
2153              if (choice==2 .and. ndgxdt==1) then
2154                jffkg=nffkge+iffkg
2155 !              Pay attention to the use of reals and imaginary parts here ...
2156 !              !$OMP CRITICAL (OPERNL4a_3)
2157                dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2158                dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg)
2159 !              !$OMP END CRITICAL (OPERNL4a_3)
2160              end if
2161              if (choice==4) then
2162                do mu=1,3
2163                  jffkg=nffkge+(iffkg-1)*9+mu
2164 !                Pay attention to the use of reals and imaginary parts here ...
2165 !                !$OMP CRITICAL (OPERNL4a_4)
2166                  dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2167                  dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
2168 !                !$OMP END CRITICAL (OPERNL4a_4)
2169                end do
2170                do mu=4,9
2171                  jffkg=nffkge+(iffkg-1)*9+mu
2172 !                Pay attention to the use of reals and imaginary parts here ...
2173 !                Also, note the multiplication by (2 pi)**2
2174 !                !$OMP CRITICAL (OPERNL4a_5)
2175                  dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg)
2176                  dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg)
2177 !                !$OMP END CRITICAL (OPERNL4a_5)
2178                end do
2179              end if
2180 !            End loop on ii=1,ilang2
2181            end do
2182 
2183            if ((choice==3.or.choice==23) .or. choice==6) then
2184 !            Compute additional tensors related to strain gradients
2185 !            ilang3 is number of unique tensor components of rank ilang+1
2186              ilang3=((ilang+2)*(ilang+3))/2
2187              jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
2188 !            Compute strain gradient tensor components
2189              do ii=1,ilang3
2190 !              Note that iffkgs is also used by ddk and 2nd derivative parts
2191                iffkgs=iffkgs+1
2192                jj=ii+jjs
2193 !              !$OMP CRITICAL (OPERNL4a_6)
2194                dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
2195                dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs)
2196 !              !$OMP END CRITICAL (OPERNL4a_6)
2197              end do
2198            end if
2199 
2200            if (choice==6) then
2201 !            Compute additional tensors related to strain 2nd derivatives
2202 !            and internal strain derivatives
2203 !            ilang6 is number of unique tensor components of rank ilang+3
2204              ilang6=((ilang+4)*(ilang+5))/2
2205              jjs=((ilang+3)*(ilang+4)*(ilang+5))/6
2206 !            Compute strain gradient tensor components
2207              do ii=1,ilang6
2208 !              Note that iffkgs is also used by ddk part
2209                iffkgs2=iffkgs2+1
2210                jj=ii+jjs
2211 !              !$OMP CRITICAL (OPERNL4a_7)
2212                d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2)
2213                d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2)
2214 !              !$OMP END CRITICAL (OPERNL4a_7)
2215              end do
2216 
2217 !            ilang4 is number of unique tensor components of rank ilang
2218              ilang4=((ilang+1)*(ilang+2))/2
2219              jjs=((ilang)*(ilang+1)*(ilang+2))/6
2220 !            Compute internal strain gradient tensor components
2221              do ii=1,ilang4
2222                iffkgs2=iffkgs2+1
2223                jj=ii+jjs
2224 !              !$OMP CRITICAL (OPERNL4a_8)
2225 !              Pay attention to the use of reals and imaginary parts here ...
2226                dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2)
2227                dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2)
2228 !              !$OMP END CRITICAL (OPERNL4a_8)
2229              end do
2230 
2231 !            ilang5 is number of unique tensor components of rank ilang+2
2232              ilang5=((ilang+3)*(ilang+4))/2
2233              jjs=((ilang+2)*(ilang+3)*(ilang+4))/6
2234 !            Compute internal strain gradient tensor components
2235              do ii=1,ilang5
2236                iffkgs2=iffkgs2+1
2237                jj=ii+jjs
2238 !              !$OMP CRITICAL (OPERNL4a_9)
2239 !              Pay attention to the use of reals and imaginary parts here ...
2240                d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2)
2241                d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2)
2242 !              !$OMP END CRITICAL (OPERNL4a_9)
2243              end do
2244            end if ! choice==6
2245 
2246            if (choice==5) then
2247 !            Compute additional tensors related to ddk with ffnl(:,2,...)
2248              ilangx=(ilang*(ilang+1))/2
2249              jjs=((ilang-1)*ilang*(ilang+1))/6
2250              do ii=1,ilangx
2251 !              Note that iffkgs is also used by strain part
2252                iffkgs=iffkgs+1
2253                jj=ii+jjs
2254 !              !$OMP CRITICAL (OPERNL4a_10)
2255                dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
2256                dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs)
2257 !              !$OMP END CRITICAL (OPERNL4a_10)
2258              end do
2259 !            Compute additional tensors related to ddk with ffnl(:,1,...)
2260              if(ilang>=2)then
2261                ilangx=((ilang-1)*ilang)/2
2262                jjs=((ilang-2)*(ilang-1)*ilang)/6
2263                do ii=1,ilangx
2264                  iffkgk=iffkgk+1
2265                  jj=ii+jjs
2266 !                !$OMP CRITICAL (OPERNL4a_11)
2267                  dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
2268                  dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk)
2269 !                !$OMP END CRITICAL (OPERNL4a_11)
2270                end do
2271              end if
2272            end if
2273 
2274 !          End projector loop
2275          end do
2276 
2277 !        End condition of non-zero projectors
2278        end if
2279 
2280 !      End angular momentum loop
2281      end do
2282 
2283 !    End loop on atoms
2284    end do
2285 
2286 !  End loop on blocks of planewaves
2287  end do
2288 !!$OMP END DO
2289 
2290  ABI_DEALLOCATE(ffkg)
2291  ABI_DEALLOCATE(kpgx)
2292  ABI_DEALLOCATE(parity)
2293  ABI_DEALLOCATE(scalars)
2294  ABI_DEALLOCATE(teffv)
2295 !!$OMP END PARALLEL
2296 
2297 !DEBUG
2298 !write(std_out,*)' opernl4a : exit'
2299 !ENDDEBUG
2300 
2301 end subroutine opernl4a

ABINIT/opernl4b [ Functions ]

[ Top ] [ Functions ]

NAME

 opernl4b

FUNCTION

 Operate with the non-local part of the hamiltonian,
 from projected quantities to reciprocal space

INPUTS

  if(choice==2 .or choice==4 .or. choice==5)
   dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
    or with respect to ddk
  if(choice==3)
   dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
  ------ Taken away in beautification because unused MS -------
   d2gxds2((2,mlang6,mincat,mproj) dummy argument here, not used
  -------------------------------------------------------------
  ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
  gmet(3,3)=metric tensor for G vecs (in bohr**-2)
  gxa(2,mlang3,mincat,mproj)= projected scalars
  ia3=gives the number of the first atom in the subset presently treated
  idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
  ispinor=1 or 2, gives the spinorial component of ffnl to be used
  ------ Taken away in beautification because unused MS -------
  istwf_k=option parameter that describes the storage of wfs
  -------------------------------------------------------------
  itypat = type of atom, needed for ffnl
  jproj(nlang)=number of projectors for each angular momentum
  kg_k(3,npw)=integer coords of planewaves in basis sphere
  kpg_k(npw,npkg)= (k+G) components and related data
  kpt(3)=real components of k point in terms of recip. translations
  lmnmax=max. number of (l,n) components over all type of psps
  matblk=dimension of the array ph3d
  mincat= maximum increment of atoms
  mlang3 = one of the dimensions of the array gxa
  mlang4 = dimension for dgxds
  ------ Taken away in beautification because unused MS -------
  mlang6 = dimension for d2gxds2
  -------------------------------------------------------------
  mproj=maximum dimension for number of projection operators for each
    angular momentum for nonlocal pseudopotential
  ndgxdt=second dimension of dgxdt
  nincat = number of atoms in the subset here treated
  nkpg=second size of array kpg_k
  nffnl=third dimension of ffnl
  nlang = number of angular momenta to be treated = 1 + highest ang. mom.
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npw  = number of plane waves in reciprocal space
  ntypat = number of type of atoms, dimension needed for ffnl
  gxa(2,mlang3,nincat,mproj)=modified projected scalars;
  NOTE that metric contractions have already been performed
  ph3d(2,npw,matblk)=three-dimensional phase factors

OUTPUT

  vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.

NOTES

 Operate with the non-local part of the hamiltonian for one type of
 atom, and within this given type of atom, for a subset of
 at most nincat atoms.

 This routine basically replaces getgla (gxa here is the former gla),
 except for the calculation of <G|dVnl/dk|C> or strain gradients.

 Present version decomposed according to iffkg
 opernl4a.f is from reciprocal space to projected quantities.

PARENTS

      nonlop_pl

CHILDREN

      dfpt_mkffkg

SOURCE

2380 subroutine opernl4b(choice,dgxds,dgxdt,ffnl,gmet,gxa,&
2381 &  ia3,idir,indlmn,ispinor,itypat,jproj,kg_k,kpg_k,kpt,&
2382 &  lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat,&
2383 &  nkpg,nlang,nloalg,npw,ntypat,ph3d,vect)
2384 
2385 
2386 !This section has been created automatically by the script Abilint (TD).
2387 !Do not modify the following lines by hand.
2388 #undef ABI_FUNC
2389 #define ABI_FUNC 'opernl4b'
2390 !End of the abilint section
2391 
2392  implicit none
2393 
2394 !Arguments ------------------------------------
2395 !scalars
2396  integer,intent(in) :: choice,ia3,idir,ispinor,itypat,lmnmax,matblk !,istwf_k
2397  integer,intent(in) :: mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat !,mlang6
2398  integer,intent(in) :: nkpg,nlang,npw,ntypat
2399 !arrays
2400  integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
2401  integer,intent(in) :: nloalg(3)
2402 !real(dp),intent(in) :: d2gxds2(2,mlang6,mincat,mproj)
2403  real(dp),intent(in) :: dgxds(2,mlang4,mincat,mproj)
2404  real(dp),intent(in) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
2405  real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3)
2406  real(dp),intent(in) :: gxa(2,mlang3,mincat,mproj),kpg_k(npw,nkpg),kpt(3)
2407  real(dp),intent(in) :: ph3d(2,npw,matblk)
2408  real(dp),intent(inout) :: vect(:,:) !vz_i
2409 
2410 !Local variables-------------------------------
2411 !scalars
2412  integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,ig,ii,ilang,ilang2,ilang3
2413  integer :: iproj,ipw,ipw1,ipw2,jffkg,jj,jump,mblkpw,nffkg
2414  integer :: nffkgd,nffkge,nffkgk,nffkgs,nincpw,nproj,ntens,start
2415  real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8
2416  real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8
2417  character(len=500) :: message
2418 !arrays
2419  integer,allocatable :: parity(:)
2420  real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
2421 
2422 ! *************************************************************************
2423 
2424 !mblkpw sets the size of blocks of planewaves
2425  mblkpw=NLO_MBLKPW
2426 
2427 !jump governs, in fine, the use of registers in the most cpu
2428 !time consuming part of the routine. Until now, jump=8 is the maximal value.
2429 !The optimal value will be machine-dependent !
2430  jump=4
2431 
2432 !Initialisation before blocking on the plane waves
2433 
2434 !Set up dimension of kpgx and allocate
2435 !ntens sets the maximum number of independent tensor components
2436 !over all allowed angular momenta; need 20 for spdf for tensors
2437 !up to rank 3; to handle stress tensor, need up to rank 5
2438  ntens=1
2439  if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5)ntens=4
2440  if(nlang>=3 .or. choice==3)ntens=10
2441  if(nlang>=4 .or. (choice==3 .and. nlang>=2) )ntens=20
2442  if(choice==3 .and. nlang>=3)ntens=35
2443  if(choice==3 .and. nlang==4)ntens=56
2444 
2445 !Set up second dimension of ffkg array, and allocate
2446  nffkg=0; nffkge=0; nffkgd=0; nffkgk=0; nffkgs=0
2447  do ilang=1,nlang
2448 !  Get the number of projectors for that angular momentum
2449    nproj=jproj(ilang)
2450 !  If there is a non-local part, accumulate the number of vectors needed
2451    if(nproj>0)then
2452      ilang2=(ilang*(ilang+1))/2
2453      nffkge=nffkge+nproj*ilang2
2454      if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
2455      if(choice==2 .or. choice==4)nffkgd=nffkgd+ndgxdt*nproj*ilang2
2456      if(choice==3)then
2457        ilang3=((ilang+2)*(ilang+3))/2
2458        nffkgs=nffkgs+nproj*ilang3
2459      end if
2460    end if
2461  end do
2462  nffkg=nffkge+nffkgd+nffkgs+nffkgk
2463 
2464 !!$OMP PARALLEL DEFAULT(PRIVATE) &
2465 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) &
2466 !!$OMP SHARED(itypat,jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang3,mlang4,mproj) &
2467 !!$OMP SHARED(ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) &
2468 !!$OMP SHARED(jump,nffkgd,nffkgk,nffkgs,mblkpw,nffkg,nffkge,ntens)
2469 
2470  ABI_ALLOCATE(ffkg,(nffkg,mblkpw))
2471  ABI_ALLOCATE(parity,(nffkg))
2472  ABI_ALLOCATE(kpgx,(mblkpw,ntens))
2473  ABI_ALLOCATE(scalars,(2,nffkg))
2474  ABI_ALLOCATE(teffv,(2,mblkpw))
2475 
2476 !Loop on subsets of plane waves (blocking)
2477 !!$OMP DO
2478  do ipw1=1,npw,mblkpw
2479 
2480    ipw2=min(npw,ipw1+mblkpw-1)
2481    nincpw=ipw2-ipw1+1
2482 
2483 !  Initialize kpgx array related to tensors defined below
2484    call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
2485 &   kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
2486 &   npw,ntens,ntypat,parity)
2487 
2488    if (choice==1 .or. choice==2 .or. choice==3 .or. choice==5) then
2489 !    Application of non-local part from projected scalars
2490 !    back to reciprocal space ...
2491 !    [this section merely computes terms which add to <G|Vnl|C>;
2492 !    nothing here is needed when various gradients are being computed]
2493 
2494 !    Loop on atoms
2495      do ia=1,nincat
2496 
2497 !      Compute the shift eventually needed to get the phases in ph3d
2498        iaph3d=ia
2499        if(nloalg(2)>0)iaph3d=ia+ia3-1
2500 
2501 !      Transfer gxa (and eventually dgxdt) in scalars with different indexing
2502        iffkg=0
2503        iffkgk=nffkge*2
2504        iffkgs=nffkge
2505        do ilang=1,nlang
2506          nproj=jproj(ilang)
2507          if (nproj>0) then
2508            ilang2=(ilang*(ilang+1))/2
2509            ilang3=((ilang+2)*(ilang+3))/2
2510            do iproj=1,nproj
2511              do ii=1,ilang2
2512                jj=ii+((ilang-1)*ilang*(ilang+1))/6
2513                iffkg=iffkg+1
2514                if(choice==1 .or. choice==3)then
2515                  scalars(1,iffkg)=gxa(1,jj,ia,iproj)
2516                  scalars(2,iffkg)=gxa(2,jj,ia,iproj)
2517                else if (choice==2 .and. ndgxdt==1) then
2518                  jffkg=nffkge+iffkg
2519 !                Pay attention to the use of reals and imaginary parts here ...
2520 !                Also, the gxa and dgxdt arrays are switched, in order
2521 !                to give the correct combination when multiplying ffkg,
2522 !                see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]]
2523                  scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
2524                  scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
2525                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
2526                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
2527                else if (choice==5) then
2528                  jffkg=nffkge+iffkg
2529 !                The gxa and dgxdt arrays are switched, in order
2530 !                to give the correct combination when multiplying ffkg,
2531                  scalars(1,jffkg)= gxa(1,jj,ia,iproj)
2532                  scalars(2,jffkg)= gxa(2,jj,ia,iproj)
2533                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
2534                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
2535                end if
2536              end do
2537              if(choice==3) then
2538                do ii=1,ilang3
2539                  iffkgs=iffkgs+1
2540                  jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
2541                  scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
2542                  scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
2543                end do
2544              end if
2545              if(ilang>=2 .and. choice==5)then
2546                do ii=1,((ilang-1)*ilang)/2
2547                  jj=ii+((ilang-2)*(ilang-1)*ilang)/6
2548                  iffkgk=iffkgk+1
2549                  scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
2550                  scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
2551                end do
2552              end if
2553            end do
2554          end if
2555        end do
2556 
2557 !      DEBUG
2558 !      write(std_out,*)' opernl4b, write scalars '
2559 !      do iffkg=1,nffkg
2560 !      write(std_out,*)iffkg,scalars(1:2,iffkg)
2561 !      end do
2562 !      ENDDEBUG
2563 
2564 !      ******* Entering the second critical part ****************************
2565 
2566 !      First, treat small nffkg; send treat the loop needed for big nffkg;
2567 !      finally treat the end of the loop needed for big nffkg
2568 
2569 !      For the time being, the maximal jump allowed is 8.
2570 
2571 !      1) Here, treat small nffkg
2572        if(nffkg<=jump)then
2573 
2574          select case(nffkg)
2575 
2576          case(1)
2577 
2578            scr1=scalars(1,1) ; sci1=scalars(2,1)
2579            ig=ipw1
2580            do ipw=1,nincpw
2581              ar=ffkg(1,ipw)*scr1 ; ai=ffkg(1,ipw)*sci1
2582              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2583              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2584              ig=ig+1
2585            end do
2586 
2587          case(2)
2588 
2589            scr1=scalars(1,1) ; sci1=scalars(2,1)
2590            scr2=scalars(1,2) ; sci2=scalars(2,2)
2591            ig=ipw1
2592            do ipw=1,nincpw
2593              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2594              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2595              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2596              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2597              ig=ig+1
2598            end do
2599 
2600          case(3)
2601 
2602            scr1=scalars(1,1) ; sci1=scalars(2,1)
2603            scr2=scalars(1,2) ; sci2=scalars(2,2)
2604            scr3=scalars(1,3) ; sci3=scalars(2,3)
2605            ig=ipw1
2606            do ipw=1,nincpw
2607              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2608              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2609              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2610              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2611              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2612              ig=ig+1
2613            end do
2614 
2615          case(4)
2616 
2617            scr1=scalars(1,1) ; sci1=scalars(2,1)
2618            scr2=scalars(1,2) ; sci2=scalars(2,2)
2619            scr3=scalars(1,3) ; sci3=scalars(2,3)
2620            scr4=scalars(1,4) ; sci4=scalars(2,4)
2621            ig=ipw1
2622            do ipw=1,nincpw
2623              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2624              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2625              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2626              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2627              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2628              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2629              ig=ig+1
2630            end do
2631 
2632          case(5)
2633 
2634            scr1=scalars(1,1) ; sci1=scalars(2,1)
2635            scr2=scalars(1,2) ; sci2=scalars(2,2)
2636            scr3=scalars(1,3) ; sci3=scalars(2,3)
2637            scr4=scalars(1,4) ; sci4=scalars(2,4)
2638            scr5=scalars(1,5) ; sci5=scalars(2,5)
2639            ig=ipw1
2640            do ipw=1,nincpw
2641              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2642              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2643              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2644              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2645              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2646              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2647              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2648              ig=ig+1
2649            end do
2650 
2651          case(6)
2652 
2653            scr1=scalars(1,1) ; sci1=scalars(2,1)
2654            scr2=scalars(1,2) ; sci2=scalars(2,2)
2655            scr3=scalars(1,3) ; sci3=scalars(2,3)
2656            scr4=scalars(1,4) ; sci4=scalars(2,4)
2657            scr5=scalars(1,5) ; sci5=scalars(2,5)
2658            scr6=scalars(1,6) ; sci6=scalars(2,6)
2659            ig=ipw1
2660            do ipw=1,nincpw
2661              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2662              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2663              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2664              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2665              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2666              ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
2667              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2668              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2669              ig=ig+1
2670            end do
2671 
2672          case(7)
2673 
2674            scr1=scalars(1,1) ; sci1=scalars(2,1)
2675            scr2=scalars(1,2) ; sci2=scalars(2,2)
2676            scr3=scalars(1,3) ; sci3=scalars(2,3)
2677            scr4=scalars(1,4) ; sci4=scalars(2,4)
2678            scr5=scalars(1,5) ; sci5=scalars(2,5)
2679            scr6=scalars(1,6) ; sci6=scalars(2,6)
2680            scr7=scalars(1,7) ; sci7=scalars(2,7)
2681            ig=ipw1
2682            do ipw=1,nincpw
2683              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2684              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2685              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2686              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2687              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2688              ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
2689              ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7
2690              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2691              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2692              ig=ig+1
2693            end do
2694 
2695          case(8)
2696 
2697            scr1=scalars(1,1) ; sci1=scalars(2,1)
2698            scr2=scalars(1,2) ; sci2=scalars(2,2)
2699            scr3=scalars(1,3) ; sci3=scalars(2,3)
2700            scr4=scalars(1,4) ; sci4=scalars(2,4)
2701            scr5=scalars(1,5) ; sci5=scalars(2,5)
2702            scr6=scalars(1,6) ; sci6=scalars(2,6)
2703            scr7=scalars(1,7) ; sci7=scalars(2,7)
2704            scr8=scalars(1,8) ; sci8=scalars(2,8)
2705            ig=ipw1
2706            do ipw=1,nincpw
2707              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
2708              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2709              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2710              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2711              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2712              ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
2713              ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7
2714              ar=ar+ffkg(8,ipw)*scr8 ; ai=ai+ffkg(8,ipw)*sci8
2715              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2716              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2717              ig=ig+1
2718            end do
2719 
2720 
2721          end select
2722 
2723        else
2724 
2725          do ipw=1,nincpw
2726            teffv(1,ipw)=0.0d0 ; teffv(2,ipw)=0.0d0
2727          end do
2728 
2729 !        2) Here treart the loop for big nffkg
2730          do start=1,nffkg-jump,jump
2731            chunk=min(jump,nffkg-jump-start+1)
2732 
2733            select case(chunk)
2734 
2735            case(1)
2736 
2737              scr1=scalars(1,start) ; sci1=scalars(2,start)
2738              do ipw=1,nincpw
2739                ar=teffv(1,ipw)            ; ai=teffv(2,ipw)
2740                ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1
2741                teffv(1,ipw)=ar            ; teffv(2,ipw)=ai
2742              end do
2743 
2744            case(2)
2745 
2746              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2747              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2748              do ipw=1,nincpw
2749                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2750                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2751                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2752                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2753              end do
2754 
2755            case(3)
2756 
2757              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2758              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2759              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2760              do ipw=1,nincpw
2761                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2762                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2763                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2764                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2765                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2766              end do
2767 
2768            case(4)
2769 
2770              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2771              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2772              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2773              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2774              do ipw=1,nincpw
2775                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2776                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2777                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2778                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2779                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2780                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2781              end do
2782 
2783            case(5)
2784 
2785              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2786              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2787              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2788              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2789              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2790              do ipw=1,nincpw
2791                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2792                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2793                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2794                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2795                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2796                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2797                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2798              end do
2799 
2800            case(6)
2801 
2802              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2803              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2804              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2805              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2806              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2807              scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2808              do ipw=1,nincpw
2809                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2810                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2811                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2812                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2813                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2814                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2815                ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2816                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2817              end do
2818 
2819            case(7)
2820 
2821              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2822              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2823              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2824              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2825              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2826              scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2827              scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2828              do ipw=1,nincpw
2829                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2830                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2831                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2832                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2833                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2834                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2835                ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2836                ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2837                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2838              end do
2839 
2840            case(8)
2841 
2842              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2843              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2844              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2845              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2846              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2847              scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2848              scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2849              scr8=scalars(1,start+7) ; sci8=scalars(2,start+7)
2850              do ipw=1,nincpw
2851                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2852                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2853                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2854                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2855                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2856                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2857                ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2858                ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2859                ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8
2860                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
2861              end do
2862 
2863            end select
2864 
2865          end do
2866 
2867 !        3) Treat the end of the loops
2868 
2869          start=nffkg-jump+1
2870 
2871          select case(jump)
2872 
2873          case(1)
2874 
2875            scr1=scalars(1,start) ; sci1=scalars(2,start)
2876            ig=ipw1
2877            do ipw=1,nincpw
2878              ar=teffv(1,ipw)            ; ai=teffv(2,ipw)
2879              ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1
2880              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2881              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2882              ig=ig+1
2883            end do
2884 
2885          case(2)
2886 
2887            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2888            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2889            ig=ipw1
2890            do ipw=1,nincpw
2891              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2892              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2893              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2894              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2895              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2896              ig=ig+1
2897            end do
2898 
2899          case(3)
2900 
2901            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2902            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2903            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2904            ig=ipw1
2905            do ipw=1,nincpw
2906              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2907              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2908              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2909              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2910              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2911              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2912              ig=ig+1
2913            end do
2914 
2915          case(4)
2916 
2917            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2918            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2919            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2920            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2921            ig=ipw1
2922            do ipw=1,nincpw
2923              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2924              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2925              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2926              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2927              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2928              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2929              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2930              ig=ig+1
2931            end do
2932 
2933          case(5)
2934 
2935            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2936            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2937            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2938            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2939            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2940            ig=ipw1
2941            do ipw=1,nincpw
2942              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2943              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2944              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2945              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2946              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2947              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2948              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2949              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2950              ig=ig+1
2951            end do
2952 
2953          case(6)
2954 
2955            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2956            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2957            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2958            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2959            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2960            scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2961            ig=ipw1
2962            do ipw=1,nincpw
2963              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2964              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2965              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2966              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2967              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2968              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2969              ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2970              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2971              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2972              ig=ig+1
2973            end do
2974 
2975          case(7)
2976 
2977            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
2978            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2979            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2980            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2981            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2982            scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2983            scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2984            ig=ipw1
2985            do ipw=1,nincpw
2986              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
2987              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
2988              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2989              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2990              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2991              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2992              ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2993              ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2994              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2995              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2996              ig=ig+1
2997            end do
2998 
2999          case(8)
3000 
3001            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
3002            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
3003            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
3004            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
3005            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
3006            scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
3007            scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
3008            scr8=scalars(1,start+7) ; sci8=scalars(2,start+7)
3009            ig=ipw1
3010            do ipw=1,nincpw
3011              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
3012              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
3013              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
3014              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
3015              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
3016              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
3017              ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
3018              ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
3019              ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8
3020              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
3021              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
3022              ig=ig+1
3023            end do
3024 
3025          end select
3026 
3027 !        End if statement for small or big nffkg
3028        end if
3029 
3030 !      ******* Leaving the critical part *********************************
3031 
3032 !      End loop on atoms
3033      end do
3034 
3035 !    End choice==1 or choice==2 or choice==3
3036    else
3037 !    Problem: choice does not make sense
3038      write(message,'(a,i0,a)' )' Input choice=',choice,' not allowed. '
3039      MSG_BUG(message)
3040    end if
3041 
3042 !  End loop on blocks of planewaves
3043  end do
3044 !!$OMP END DO
3045 
3046  ABI_DEALLOCATE(ffkg)
3047  ABI_DEALLOCATE(kpgx)
3048  ABI_DEALLOCATE(parity)
3049  ABI_DEALLOCATE(scalars)
3050  ABI_DEALLOCATE(teffv)
3051 !!$OMP END PARALLEL
3052 
3053 end subroutine opernl4b