TABLE OF CONTENTS


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)

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, DRH)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

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

OUTPUT

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

NOTES

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

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

PARENTS

      nonlop_pl

CHILDREN

      dfpt_mkffkg

SOURCE

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