TABLE OF CONTENTS


ABINIT/m_opernl [ Modules ]

[ Top ] [ Modules ]

NAME

  m_opernl

FUNCTION

COPYRIGHT

 Copyright (C) 1998-2022 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_opernl
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27 
28  use m_mkffkg, only : mkffkg, dfpt_mkffkg
29 
30  implicit none
31 
32  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.

SOURCE

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

SOURCE

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

SOURCE

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

SOURCE

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