TABLE OF CONTENTS


ABINIT/opernl4b [ Functions ]

[ Top ] [ Functions ]

NAME

 opernl4b

FUNCTION

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

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

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

 85 #if defined HAVE_CONFIG_H
 86 #include "config.h"
 87 #endif
 88 
 89 #include "abi_common.h"
 90 
 91 subroutine opernl4b(choice,dgxds,dgxdt,ffnl,gmet,gxa,&
 92 &  ia3,idir,indlmn,ispinor,itypat,jproj,kg_k,kpg_k,kpt,&
 93 &  lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat,&
 94 &  nkpg,nlang,nloalg,npw,ntypat,ph3d,vect)
 95 
 96  use defs_basis
 97  use m_profiling_abi
 98  use m_errors
 99 
100 !This section has been created automatically by the script Abilint (TD).
101 !Do not modify the following lines by hand.
102 #undef ABI_FUNC
103 #define ABI_FUNC 'opernl4b'
104  use interfaces_66_nonlocal, except_this_one => opernl4b
105 !End of the abilint section
106 
107  implicit none
108 
109 !Arguments ------------------------------------
110 !scalars
111  integer,intent(in) :: choice,ia3,idir,ispinor,itypat,lmnmax,matblk !,istwf_k
112  integer,intent(in) :: mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat !,mlang6
113  integer,intent(in) :: nkpg,nlang,npw,ntypat
114 !arrays
115  integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
116  integer,intent(in) :: nloalg(3)
117 !real(dp),intent(in) :: d2gxds2(2,mlang6,mincat,mproj)
118  real(dp),intent(in) :: dgxds(2,mlang4,mincat,mproj)
119  real(dp),intent(in) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
120  real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3)
121  real(dp),intent(in) :: gxa(2,mlang3,mincat,mproj),kpg_k(npw,nkpg),kpt(3)
122  real(dp),intent(in) :: ph3d(2,npw,matblk)
123  real(dp),intent(inout) :: vect(:,:) !vz_i
124 
125 !Local variables-------------------------------
126 !scalars
127  integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,ig,ii,ilang,ilang2,ilang3
128  integer :: iproj,ipw,ipw1,ipw2,jffkg,jj,jump,mblkpw,nffkg
129  integer :: nffkgd,nffkge,nffkgk,nffkgs,nincpw,nproj,ntens,start
130  real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8
131  real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8
132  character(len=500) :: message
133 !arrays
134  integer,allocatable :: parity(:)
135  real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
136 
137 ! *************************************************************************
138 
139 !mblkpw sets the size of blocks of planewaves
140  mblkpw=NLO_MBLKPW
141 
142 !jump governs, in fine, the use of registers in the most cpu
143 !time consuming part of the routine. Until now, jump=8 is the maximal value.
144 !The optimal value will be machine-dependent !
145  jump=4
146 
147 !Initialisation before blocking on the plane waves
148 
149 !Set up dimension of kpgx and allocate
150 !ntens sets the maximum number of independent tensor components
151 !over all allowed angular momenta; need 20 for spdf for tensors
152 !up to rank 3; to handle stress tensor, need up to rank 5
153  ntens=1
154  if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5)ntens=4
155  if(nlang>=3 .or. choice==3)ntens=10
156  if(nlang>=4 .or. (choice==3 .and. nlang>=2) )ntens=20
157  if(choice==3 .and. nlang>=3)ntens=35
158  if(choice==3 .and. nlang==4)ntens=56
159 
160 !Set up second dimension of ffkg array, and allocate
161  nffkg=0; nffkge=0; nffkgd=0; nffkgk=0; nffkgs=0
162  do ilang=1,nlang
163 !  Get the number of projectors for that angular momentum
164    nproj=jproj(ilang)
165 !  If there is a non-local part, accumulate the number of vectors needed
166    if(nproj>0)then
167      ilang2=(ilang*(ilang+1))/2
168      nffkge=nffkge+nproj*ilang2
169      if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
170      if(choice==2 .or. choice==4)nffkgd=nffkgd+ndgxdt*nproj*ilang2
171      if(choice==3)then
172        ilang3=((ilang+2)*(ilang+3))/2
173        nffkgs=nffkgs+nproj*ilang3
174      end if
175    end if
176  end do
177  nffkg=nffkge+nffkgd+nffkgs+nffkgk
178 
179 !!$OMP PARALLEL DEFAULT(PRIVATE) &
180 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) &
181 !!$OMP SHARED(itypat,jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang3,mlang4,mproj) &
182 !!$OMP SHARED(ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) &
183 !!$OMP SHARED(jump,nffkgd,nffkgk,nffkgs,mblkpw,nffkg,nffkge,ntens)
184 
185  ABI_ALLOCATE(ffkg,(nffkg,mblkpw))
186  ABI_ALLOCATE(parity,(nffkg))
187  ABI_ALLOCATE(kpgx,(mblkpw,ntens))
188  ABI_ALLOCATE(scalars,(2,nffkg))
189  ABI_ALLOCATE(teffv,(2,mblkpw))
190 
191 !Loop on subsets of plane waves (blocking)
192 !!$OMP DO 
193  do ipw1=1,npw,mblkpw
194 
195    ipw2=min(npw,ipw1+mblkpw-1)
196    nincpw=ipw2-ipw1+1
197 
198 !  Initialize kpgx array related to tensors defined below
199    call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
200 &   kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
201 &   npw,ntens,ntypat,parity)
202 
203    if (choice==1 .or. choice==2 .or. choice==3 .or. choice==5) then
204 !    Application of non-local part from projected scalars
205 !    back to reciprocal space ...
206 !    [this section merely computes terms which add to <G|Vnl|C>;
207 !    nothing here is needed when various gradients are being computed]
208 
209 !    Loop on atoms
210      do ia=1,nincat
211 
212 !      Compute the shift eventually needed to get the phases in ph3d
213        iaph3d=ia
214        if(nloalg(2)>0)iaph3d=ia+ia3-1
215 
216 !      Transfer gxa (and eventually dgxdt) in scalars with different indexing
217        iffkg=0
218        iffkgk=nffkge*2
219        iffkgs=nffkge
220        do ilang=1,nlang
221          nproj=jproj(ilang)
222          if (nproj>0) then
223            ilang2=(ilang*(ilang+1))/2
224            ilang3=((ilang+2)*(ilang+3))/2
225            do iproj=1,nproj
226              do ii=1,ilang2
227                jj=ii+((ilang-1)*ilang*(ilang+1))/6
228                iffkg=iffkg+1
229                if(choice==1 .or. choice==3)then
230                  scalars(1,iffkg)=gxa(1,jj,ia,iproj)
231                  scalars(2,iffkg)=gxa(2,jj,ia,iproj)
232                else if (choice==2 .and. ndgxdt==1) then
233                  jffkg=nffkge+iffkg
234 !                Pay attention to the use of reals and imaginary parts here ...
235 !                Also, the gxa and dgxdt arrays are switched, in order
236 !                to give the correct combination when multiplying ffkg,
237 !                see Eq.(53) of PRB55,10337(1997)
238                  scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
239                  scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
240                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
241                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
242                else if (choice==5) then
243                  jffkg=nffkge+iffkg
244 !                The gxa and dgxdt arrays are switched, in order
245 !                to give the correct combination when multiplying ffkg,
246                  scalars(1,jffkg)= gxa(1,jj,ia,iproj)
247                  scalars(2,jffkg)= gxa(2,jj,ia,iproj)
248                  scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
249                  scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
250                end if
251              end do
252              if(choice==3) then
253                do ii=1,ilang3
254                  iffkgs=iffkgs+1
255                  jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
256                  scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
257                  scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
258                end do
259              end if
260              if(ilang>=2 .and. choice==5)then
261                do ii=1,((ilang-1)*ilang)/2
262                  jj=ii+((ilang-2)*(ilang-1)*ilang)/6
263                  iffkgk=iffkgk+1
264                  scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
265                  scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
266                end do
267              end if
268            end do
269          end if
270        end do
271 
272 !      DEBUG
273 !      write(std_out,*)' opernl4b, write scalars '
274 !      do iffkg=1,nffkg
275 !      write(std_out,*)iffkg,scalars(1:2,iffkg)
276 !      end do
277 !      ENDDEBUG
278 
279 !      ******* Entering the second critical part ****************************
280 
281 !      First, treat small nffkg; send treat the loop needed for big nffkg;
282 !      finally treat the end of the loop needed for big nffkg
283 
284 !      For the time being, the maximal jump allowed is 8.
285 
286 !      1) Here, treat small nffkg
287        if(nffkg<=jump)then
288 
289          select case(nffkg)
290 
291          case(1)
292 
293            scr1=scalars(1,1) ; sci1=scalars(2,1)
294            ig=ipw1
295            do ipw=1,nincpw
296              ar=ffkg(1,ipw)*scr1 ; ai=ffkg(1,ipw)*sci1
297              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
298              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
299              ig=ig+1
300            end do
301 
302          case(2)
303 
304            scr1=scalars(1,1) ; sci1=scalars(2,1)
305            scr2=scalars(1,2) ; sci2=scalars(2,2)
306            ig=ipw1
307            do ipw=1,nincpw
308              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
309              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
310              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
311              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
312              ig=ig+1
313            end do
314 
315          case(3)
316 
317            scr1=scalars(1,1) ; sci1=scalars(2,1)
318            scr2=scalars(1,2) ; sci2=scalars(2,2)
319            scr3=scalars(1,3) ; sci3=scalars(2,3)
320            ig=ipw1
321            do ipw=1,nincpw
322              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
323              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
324              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
325              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
326              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
327              ig=ig+1
328            end do
329 
330          case(4)
331 
332            scr1=scalars(1,1) ; sci1=scalars(2,1)
333            scr2=scalars(1,2) ; sci2=scalars(2,2)
334            scr3=scalars(1,3) ; sci3=scalars(2,3)
335            scr4=scalars(1,4) ; sci4=scalars(2,4)
336            ig=ipw1
337            do ipw=1,nincpw
338              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
339              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
340              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
341              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
342              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
343              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
344              ig=ig+1
345            end do
346 
347          case(5)
348 
349            scr1=scalars(1,1) ; sci1=scalars(2,1)
350            scr2=scalars(1,2) ; sci2=scalars(2,2)
351            scr3=scalars(1,3) ; sci3=scalars(2,3)
352            scr4=scalars(1,4) ; sci4=scalars(2,4)
353            scr5=scalars(1,5) ; sci5=scalars(2,5)
354            ig=ipw1
355            do ipw=1,nincpw
356              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
357              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
358              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
359              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
360              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
361              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
362              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
363              ig=ig+1
364            end do
365 
366          case(6)
367 
368            scr1=scalars(1,1) ; sci1=scalars(2,1)
369            scr2=scalars(1,2) ; sci2=scalars(2,2)
370            scr3=scalars(1,3) ; sci3=scalars(2,3)
371            scr4=scalars(1,4) ; sci4=scalars(2,4)
372            scr5=scalars(1,5) ; sci5=scalars(2,5)
373            scr6=scalars(1,6) ; sci6=scalars(2,6)
374            ig=ipw1
375            do ipw=1,nincpw
376              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
377              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
378              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
379              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
380              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
381              ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
382              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
383              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
384              ig=ig+1
385            end do
386 
387          case(7)
388 
389            scr1=scalars(1,1) ; sci1=scalars(2,1)
390            scr2=scalars(1,2) ; sci2=scalars(2,2)
391            scr3=scalars(1,3) ; sci3=scalars(2,3)
392            scr4=scalars(1,4) ; sci4=scalars(2,4)
393            scr5=scalars(1,5) ; sci5=scalars(2,5)
394            scr6=scalars(1,6) ; sci6=scalars(2,6)
395            scr7=scalars(1,7) ; sci7=scalars(2,7)
396            ig=ipw1
397            do ipw=1,nincpw
398              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
399              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
400              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
401              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
402              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
403              ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
404              ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7
405              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
406              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
407              ig=ig+1
408            end do
409 
410          case(8)
411 
412            scr1=scalars(1,1) ; sci1=scalars(2,1)
413            scr2=scalars(1,2) ; sci2=scalars(2,2)
414            scr3=scalars(1,3) ; sci3=scalars(2,3)
415            scr4=scalars(1,4) ; sci4=scalars(2,4)
416            scr5=scalars(1,5) ; sci5=scalars(2,5)
417            scr6=scalars(1,6) ; sci6=scalars(2,6)
418            scr7=scalars(1,7) ; sci7=scalars(2,7)
419            scr8=scalars(1,8) ; sci8=scalars(2,8)
420            ig=ipw1
421            do ipw=1,nincpw
422              ar=   ffkg(1,ipw)*scr1 ; ai=   ffkg(1,ipw)*sci1
423              ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
424              ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
425              ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
426              ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
427              ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
428              ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7
429              ar=ar+ffkg(8,ipw)*scr8 ; ai=ai+ffkg(8,ipw)*sci8
430              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
431              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
432              ig=ig+1
433            end do
434 
435 
436          end select
437 
438        else
439 
440          do ipw=1,nincpw
441            teffv(1,ipw)=0.0d0 ; teffv(2,ipw)=0.0d0
442          end do
443 
444 !        2) Here treart the loop for big nffkg
445          do start=1,nffkg-jump,jump
446            chunk=min(jump,nffkg-jump-start+1)
447 
448            select case(chunk)
449 
450            case(1)
451 
452              scr1=scalars(1,start) ; sci1=scalars(2,start)
453              do ipw=1,nincpw
454                ar=teffv(1,ipw)            ; ai=teffv(2,ipw)
455                ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1
456                teffv(1,ipw)=ar            ; teffv(2,ipw)=ai
457              end do
458 
459            case(2)
460 
461              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
462              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
463              do ipw=1,nincpw
464                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
465                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
466                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
467                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
468              end do
469 
470            case(3)
471 
472              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
473              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
474              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
475              do ipw=1,nincpw
476                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
477                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
478                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
479                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
480                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
481              end do
482 
483            case(4)
484 
485              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
486              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
487              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
488              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
489              do ipw=1,nincpw
490                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
491                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
492                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
493                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
494                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
495                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
496              end do
497 
498            case(5)
499 
500              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
501              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
502              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
503              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
504              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
505              do ipw=1,nincpw
506                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
507                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
508                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
509                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
510                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
511                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
512                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
513              end do
514 
515            case(6)
516 
517              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
518              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
519              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
520              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
521              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
522              scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
523              do ipw=1,nincpw
524                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
525                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
526                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
527                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
528                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
529                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
530                ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
531                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
532              end do
533 
534            case(7)
535 
536              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
537              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
538              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
539              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
540              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
541              scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
542              scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
543              do ipw=1,nincpw
544                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
545                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
546                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
547                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
548                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
549                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
550                ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
551                ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
552                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
553              end do
554 
555            case(8)
556 
557              scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
558              scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
559              scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
560              scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
561              scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
562              scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
563              scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
564              scr8=scalars(1,start+7) ; sci8=scalars(2,start+7)
565              do ipw=1,nincpw
566                ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
567                ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
568                ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
569                ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
570                ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
571                ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
572                ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
573                ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
574                ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8
575                teffv(1,ipw)=ar              ; teffv(2,ipw)=ai
576              end do
577 
578            end select
579 
580          end do
581 
582 !        3) Treat the end of the loops
583 
584          start=nffkg-jump+1
585 
586          select case(jump)
587 
588          case(1)
589 
590            scr1=scalars(1,start) ; sci1=scalars(2,start)
591            ig=ipw1
592            do ipw=1,nincpw
593              ar=teffv(1,ipw)            ; ai=teffv(2,ipw)
594              ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1
595              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
596              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
597              ig=ig+1
598            end do
599 
600          case(2)
601 
602            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
603            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
604            ig=ipw1
605            do ipw=1,nincpw
606              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
607              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
608              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
609              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
610              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
611              ig=ig+1
612            end do
613 
614          case(3)
615 
616            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
617            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
618            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
619            ig=ipw1
620            do ipw=1,nincpw
621              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
622              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
623              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
624              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
625              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
626              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
627              ig=ig+1
628            end do
629 
630          case(4)
631 
632            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
633            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
634            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
635            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
636            ig=ipw1
637            do ipw=1,nincpw
638              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
639              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
640              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
641              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
642              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
643              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
644              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
645              ig=ig+1
646            end do
647 
648          case(5)
649 
650            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
651            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
652            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
653            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
654            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
655            ig=ipw1
656            do ipw=1,nincpw
657              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
658              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
659              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
660              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
661              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
662              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
663              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
664              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
665              ig=ig+1
666            end do
667 
668          case(6)
669 
670            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
671            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
672            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
673            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
674            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
675            scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
676            ig=ipw1
677            do ipw=1,nincpw
678              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
679              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
680              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
681              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
682              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
683              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
684              ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
685              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
686              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
687              ig=ig+1
688            end do
689 
690          case(7)
691 
692            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
693            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
694            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
695            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
696            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
697            scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
698            scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
699            ig=ipw1
700            do ipw=1,nincpw
701              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
702              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
703              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
704              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
705              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
706              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
707              ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
708              ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
709              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
710              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
711              ig=ig+1
712            end do
713 
714          case(8)
715 
716            scr1=scalars(1,start  ) ; sci1=scalars(2,start  )
717            scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
718            scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
719            scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
720            scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
721            scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
722            scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
723            scr8=scalars(1,start+7) ; sci8=scalars(2,start+7)
724            ig=ipw1
725            do ipw=1,nincpw
726              ar=teffv(1,ipw)              ; ai=teffv(2,ipw)
727              ar=ar+ffkg(start  ,ipw)*scr1 ; ai=ai+ffkg(start  ,ipw)*sci1
728              ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
729              ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
730              ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
731              ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
732              ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
733              ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
734              ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8
735              vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
736              vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
737              ig=ig+1
738            end do
739 
740          end select
741 
742 !        End if statement for small or big nffkg
743        end if
744 
745 !      ******* Leaving the critical part *********************************
746 
747 !      End loop on atoms
748      end do
749 
750 !    End choice==1 or choice==2 or choice==3
751    else
752 !    Problem: choice does not make sense
753      write(message,'(a,i0,a)' )' Input choice=',choice,' not allowed. '
754      MSG_BUG(message)
755    end if
756 
757 !  End loop on blocks of planewaves
758  end do
759 !!$OMP END DO
760 
761  ABI_DEALLOCATE(ffkg)
762  ABI_DEALLOCATE(kpgx)
763  ABI_DEALLOCATE(parity)
764  ABI_DEALLOCATE(scalars)
765  ABI_DEALLOCATE(teffv)
766 !!$OMP END PARALLEL
767 
768 end subroutine opernl4b