TABLE OF CONTENTS


ABINIT/m_opernld_ylm [ Modules ]

[ Top ] [ Modules ]

NAME

  m_opernld_ylm

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (MT)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_opernld_ylm
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  implicit none
34 
35  private

ABINIT/opernld_ylm [ Functions ]

[ Top ] [ Functions ]

NAME

 opernld_ylm

FUNCTION

 * Operate with the non-local part of the hamiltonian,
   in order to get contributions to energy/forces/stress/dyn.matrix/elst tens.
   from projected scalars
 * Operate with the non-local projectors and the overlap matrix Sij
   in order to get contributions to <c|S|c>
   from projected scalars

INPUTS

  choice=chooses possible output
  cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
        2 if <p_lmn|c> scalars are complex
  cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
  d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)=2nd gradients of projected scalars
  dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)=gradients of projected scalars
  dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
                                                    related to Vnl (NL operator)
  dgxdtfac_sij(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
                                                        related to Sij (overlap)
  gx(cplex,nlmn,nincat,nspinor)= projected scalars
  gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
  gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
  ia3=gives the absolute number of the first atom in the subset presently treated
  natom=number of atoms in cell
  nd2gxdt=second dimension of d2gxdt
  ndgxdt=second dimension of dgxdt
  ndgxdtfac=second dimension of dgxdtfac
  nincat=number of atoms in the subset here treated
  nlmn=number of (l,m,n) numbers for current type of atom
  nnlout=dimension of enlout
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  paw_opt= define the nonlocal operator concerned with:
           paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
           paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
           paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
           paw_opt=3 : PAW overlap matrix (Sij)
           paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)

OUTPUT

  (see side effects)

SIDE EFFECTS

 --If (paw_opt==0, 1 or 2)
    enlout(nnlout)= contribution to the non-local part of the following properties:
      if choice=1 : enlout(1)             -> the energy
      if choice=2 : enlout(3*natom)       -> 1st deriv. of energy wrt atm. pos (forces)
      if choice=3 : enlout(6)             -> 1st deriv. of energy wrt strain (stresses)
      if choice=4 : enlout(6*natom)       -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
      if choice=23: enlout(6+3*natom)     -> 1st deriv. of energy wrt atm. pos (forces) and
                                             1st deriv. of energy wrt strain (stresses)
      if choice=24: enlout(9*natom)       -> 1st deriv. of energy wrt atm. pos (forces) and
                                             2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
      if choice=5 : enlout(3)             -> 1st deriv. of energy wrt k
      if choice=53: enlout(3)             -> 1st deriv. (twist) of energy wrt k
      if choice=54: enlout(18*natom)      -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
      if choice=55: enlout(36)            -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
      if choice=6 : enlout(36+18*natom)   -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
                                             2nd deriv. of energy wrt to atm. pos and strain (internal strain)
      if choice=8 : enlout(6)             -> 2nd deriv. of energy wrt 2 k
      if choice=81: enlout(18)            -> 2nd deriv. of energy wrt k and right k
 --If (paw_opt==3)
      if choice=1 : enlout(1)             -> contribution to <c|S|c> (note: not including <c|c>)
      if choice=2 : enlout(3*natom)       -> contribution to <c|dS/d_atm.pos|c>
      if choice=54: enlout(18*natom)      -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
      if choice=55: enlout(36)            -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
      if choice=8 : enlout(6)             -> 2nd deriv. of energy wrt 2 k
      if choice=81: enlout(18)            -> 2nd deriv. of energy wrt k and right k
 --If (paw_opt==4)
      not available

NOTES

 Operate for one type of atom, and within this given type of atom,
 for a subset of at most nincat atoms.

PARENTS

      nonlop_ylm

CHILDREN

SOURCE

129 subroutine opernld_ylm(choice,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
130 &                      enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,nd2gxdt,ndgxdt,&
131 &                      ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt,strnlk)
132 
133 
134 !This section has been created automatically by the script Abilint (TD).
135 !Do not modify the following lines by hand.
136 #undef ABI_FUNC
137 #define ABI_FUNC 'opernld_ylm'
138 !End of the abilint section
139 
140  implicit none
141 
142 !Arguments ------------------------------------
143 !scalars
144  integer,intent(in) :: choice,cplex,cplex_fac,ia3,natom,nd2gxdt,ndgxdt
145  integer,intent(in) :: ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt
146  real(dp),intent(inout) :: enlk
147 !arrays
148  real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)
149  real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)
150  real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)
151  real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3))
152  real(dp),intent(in) :: gx(cplex,nlmn,nincat,nspinor),gxfac(cplex_fac,nlmn,nincat,nspinor)
153  real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
154  real(dp),intent(inout) :: ddkk(6),enlout(nnlout),fnlk(3*natom),strnlk(6)
155 
156 !Local variables-------------------------------
157 !scalars
158  integer :: ia,iashift,ilmn,iplex,ishift,ispinor,mu,mua,mua1,mua2,mub,mushift,mut,muu,nu,nushift
159  real(dp) :: dummy
160 !arrays
161  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
162  integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
163  integer,parameter :: twist_dir(6)=(/2,3,3,1,1,2/)
164  real(dp) :: d2gx(cplex),enlj(6*cplex),gxfacj(cplex)
165  real(dp),allocatable :: enljj(:)
166  complex(dpc),allocatable :: cft(:,:), cfu(:,:)
167 
168 ! *************************************************************************
169 
170  ABI_CHECK(cplex_fac>=cplex,'BUG: invalid cplex_fac<cplex!')
171 
172  if (paw_opt==0.or.paw_opt==1.or.paw_opt==2) then
173 
174 !  ============== Accumulate the non-local energy ===============
175    if (choice==1) then
176      do ispinor=1,nspinor
177        do ia=1,nincat
178          do ilmn=1,nlmn
179            do iplex=1,cplex
180              enlout(1)=enlout(1)+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor)
181            end do
182          end do
183        end do
184      end do
185    end if
186 
187 !  ============ Accumulate the forces contributions =============
188    if (choice==2.or.choice==23.or.choice==24) then
189      ishift=0;if (choice==23) ishift=6
190      do ispinor=1,nspinor
191        do ia=1,nincat
192          enlj(1:3)=zero
193          iashift=3*(ia+ia3-2)+ishift
194          do ilmn=1,nlmn
195            do mu=1,3
196              dummy = zero  ! Dummy needed here to get the correct forces with intel -O3
197              do iplex=1,cplex
198                !enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
199                dummy=dummy+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
200              end do
201              enlj(mu)=enlj(mu)+dummy
202            end do
203          end do
204          enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3)
205        end do
206      end do
207    end if
208 
209 !  ======== Accumulate the stress tensor contributions ==========
210    if (choice==3.or.choice==23) then
211      enlj(1:6)=zero
212      do ispinor=1,nspinor
213        do ia=1,nincat
214          do ilmn=1,nlmn
215            gxfacj(1:cplex)=gxfac(1:cplex,ilmn,ia,ispinor)
216            do iplex=1,cplex
217              enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor)
218            end do
219            do mu=1,6
220              do iplex=1,cplex
221                enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor)
222              end do
223            end do
224          end do
225        end do
226      end do
227      enlout(1:6)=enlout(1:6)+two*enlj(1:6)
228    end if
229 
230 !  ====== Accumulate the dynamical matrix contributions =========
231    if (choice==4.or.choice==24) then
232      ishift=0;if (choice==24) ishift=3*natom
233      do ispinor=1,nspinor
234        do ia=1,nincat
235          enlj(1:6)=zero
236          iashift=6*(ia+ia3-2)+ishift
237          do ilmn=1,nlmn
238            do mu=1,6
239              mua=alpha(mu);mub=beta(mu)
240              do iplex=1,cplex
241                enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
242 &               +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,mua,ilmn,ia,ispinor)
243              end do
244            end do
245          end do
246          enlout(iashift+1:iashift+6)=enlout(iashift+1:iashift+6)+two*enlj(1:6)
247        end do
248      end do
249    end if
250 
251 !  ======== Accumulate the contributions of derivatives of E wrt to k ==========
252    if (choice==5) then
253      enlj(1:3)=zero
254      do ispinor=1,nspinor
255        do ia=1,nincat
256          if(cplex==2)then
257            do ilmn=1,nlmn
258              do mu=1,3
259                do iplex=1,cplex
260                  enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor)
261                end do
262              end do
263            end do
264 !        If cplex=1, dgxdt is pure imaginary; thus there is no contribution
265          else if (cplex_fac==2) then
266            do ilmn=1,nlmn
267              do mu=1,3
268                enlj(mu)=enlj(mu)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
269              end do
270            end do
271          end if
272        end do
273      end do
274      enlout(1:3)=enlout(1:3)+two*enlj(1:3)
275    end if
276 
277 !  ======== Accumulate the contributions of partial derivatives of E wrt to k ==========
278 !  Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k
279    if (choice==51.or.choice==52) then
280      enlj(1:6)=zero
281      do ispinor=1,nspinor
282        do ia=1,nincat
283          if(cplex==2)then
284            do ilmn=1,nlmn
285              do mu=1,3
286                enlj(2*mu-1)=enlj(2*mu-1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) &
287 &               +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor)
288                enlj(2*mu  )=enlj(2*mu  )+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) &
289 &               -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
290              end do
291            end do
292          else if (cplex_fac==2) then
293            do ilmn=1,nlmn
294              do mu=1,3
295                enlj(2*mu-1)=enlj(2*mu-1)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
296                enlj(2*mu  )=enlj(2*mu  )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
297              end do
298            end do
299          else if (cplex_fac==1) then
300            do ilmn=1,nlmn
301              do mu=1,3
302                enlj(2*mu  )=enlj(2*mu  )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
303              end do
304            end do
305          end if
306        end do
307      end do
308      if (choice==52) then
309        enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6)
310      end if
311      enlout(1:6)=enlout(1:6)+enlj(1:6)
312    end if
313 
314 !  ======== Accumulate the contributions of twist derivatives of E wrt to k ==========
315    if (choice==53) then
316      enlj(1:3)=zero
317      ABI_ALLOCATE(cft,(3,nlmn))
318      ABI_ALLOCATE(cfu,(3,nlmn))
319 !    If cplex=1, dgxdt is pure imaginary;
320 !    If cplex_fac=1, dgxdtfac is pure imaginary;
321      do ispinor=1,nspinor
322        do ia=1,nincat
323          if(cplex==2)then
324            cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
325          else
326            cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
327          end if
328          if(cplex_fac==2)then
329            cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor))
330          else
331            cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor))
332          end if
333          do ilmn=1,nlmn
334            do mu=1,3
335              mut = twist_dir(2*mu-1)
336              muu = twist_dir(2*mu)
337              enlj(mu) = enlj(mu) + aimag(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
338              enlj(mu) = enlj(mu) - aimag(conjg(cft(muu,ilmn))*cfu(mut,ilmn))
339            end do
340          end do
341        end do
342      end do
343      enlout(1:3)=enlout(1:3)+enlj(1:3)
344      ABI_DEALLOCATE(cft)
345      ABI_DEALLOCATE(cfu)
346    end if
347 
348 !  ====== Accumulate the effective charges contributions =========
349    if (choice==54) then
350      ABI_ALLOCATE(enljj,(18))
351      do ispinor=1,nspinor
352        do ia=1,nincat
353          enljj(1:18)=zero
354          iashift=18*(ia+ia3-2)
355 !        If cplex=1, dgxdt is real for atm. pos, pure imaginary for k;
356 !        If cplex_fac=1, dgxdtfac is pure imaginary for k;
357          if(cplex==2.and.cplex_fac==2) then
358            do ilmn=1,nlmn
359              mu=1;nu=1
360              do mua=1,3 ! atm. pos
361                do mub=1,3 ! k
362                  enljj(nu)=enljj(nu) &
363 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
364 &                 +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) &
365 &                 +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) &
366 &                 +gxfac(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor)
367                  enljj(nu+1)=enljj(nu+1) &
368 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) &
369 &                 -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
370 &                 +gxfac(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) &
371 &                 -gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
372                  mu=mu+1;nu=nu+2
373                end do
374              end do
375            end do
376          else if(cplex==1.and.cplex_fac==2)then
377            do ilmn=1,nlmn
378              mu=1;nu=1
379              do mua=1,3 ! atm. pos
380                do mub=1,3 ! k
381                  enljj(nu)=enljj(nu) &
382 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
383 &                 +gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
384                  enljj(nu+1)=enljj(nu+1) &
385 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) &
386 &                 +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
387                  mu=mu+1;nu=nu+2
388                end do
389              end do
390            end do
391          else if(cplex==1.and.cplex_fac==1)then
392            do ilmn=1,nlmn
393              mu=1;nu=1
394              do mua=1,3 ! atm. pos
395                do mub=1,3 ! k
396                  enljj(nu+1)=enljj(nu+1) &
397 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
398 &                 +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
399                  mu=mu+1;nu=nu+2
400                end do
401              end do
402            end do
403          end if
404          enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18)
405        end do
406      end do
407      ABI_DEALLOCATE(enljj)
408    end if
409 
410 !  ====== Accumulate the piezoelectric tensor contributions =========
411    if (choice==55) then
412      ABI_ALLOCATE(enljj,(36))
413      do ispinor=1,nspinor
414        do ia=1,nincat
415          enljj(1:36)=zero;enlj(:)=zero
416 !        If cplex=1, dgxdt is real for strain, pure imaginary for k;
417 !        If cplex_fac=1, dgxdtfac is pure imaginary for k;
418          if(cplex==2.and.cplex_fac==2) then
419            do ilmn=1,nlmn
420 !            First compute 2nd-derivative contribution
421              mu=1
422              do mua=1,6 ! strain (lambda,nu)
423                mua1=alpha(mua) ! (nu)
424                mua2=beta(mua)  ! (lambda)
425                do mub=1,3 ! k (mu)
426                  muu=3*(gamma(mua1,mub)-1)+mua2
427                  mut=3*(gamma(mua2,mub)-1)+mua1
428                  d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) &
429 &                 +d2gxdt(1:cplex,mut,ilmn,ia,ispinor))
430                  enljj(mu)=enljj(mu) &
431 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
432 &                 +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) &
433 &                 +gxfac(1,ilmn,ia,ispinor)*d2gx(1)+gxfac(2,ilmn,ia,ispinor)*d2gx(2)
434                  enljj(mu+1)=enljj(mu+1) &
435 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) &
436 &                 -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
437 &                 +gxfac(1,ilmn,ia,ispinor)*d2gx(2)-gxfac(2,ilmn,ia,ispinor)*d2gx(1)
438                  mu=mu+2
439                end do
440              end do
441 !            Then store 1st-derivative contribution
442              mu=1
443              do nu=1,3
444                enlj(mu  )=enlj(mu  )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) &
445 &               +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor)
446                enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) &
447 &               -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
448                mu=mu+2
449              end do
450            end do
451          else if(cplex==1.and.cplex_fac==2)then
452            do ilmn=1,nlmn
453 !            First compute 2nd-derivative contribution
454              mu=1
455              do mua=1,6 ! strain (lambda,nu)
456                mua1=alpha(mua) ! (nu)
457                mua2=beta(mua)  ! (lambda)
458                do mub=1,3 ! k (mu)
459                  muu=3*(gamma(mua1,mub)-1)+mua2
460                  mut=3*(gamma(mua2,mub)-1)+mua1
461                  d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor))
462                  enljj(mu)=enljj(mu) &
463 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
464 &                 +gxfac(2,ilmn,ia,ispinor)*d2gx(1)
465                  enljj(mu+1)=enljj(mu+1) &
466 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) &
467 &                 +gxfac(1,ilmn,ia,ispinor)*d2gx(1)
468                  mu=mu+2
469                end do
470              end do
471 !            Then store 1st-derivative contribution
472              mu=1
473              do nu=1,3
474                enlj(mu  )=enlj(mu  )+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
475                enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
476                mu=mu+2
477              end do
478            end do
479          else if(cplex==1.and.cplex_fac==1)then
480            do ilmn=1,nlmn
481              mu=1
482              do mua=1,6 ! strain (lambda,nu)
483                mua1=alpha(mua) ! (nu)
484                mua2=beta(mua)  ! (lambda)
485                do mub=1,3 ! k (mu)
486                  muu=3*(gamma(mua1,mub)-1)+mua2
487                  mut=3*(gamma(mua2,mub)-1)+mua1
488                  d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor))
489                  enljj(mu+1)=enljj(mu+1) &
490 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
491 &                 +gxfac(1,ilmn,ia,ispinor)*d2gx(1)
492                  mu=mu+2
493                end do
494              end do
495 !            Then store 1st-derivative contribution
496              mu=1
497              do nu=1,3
498                enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
499                mu=mu+2
500              end do
501            end do
502          end if
503          enlout(1:36)=enlout(1:36)+enljj(1:36)
504          ddkk(1:6)=ddkk(1:6)+enlj(1:6)
505        end do
506      end do
507      ABI_DEALLOCATE(enljj)
508    end if
509 
510 !  ======= Accumulate the elastic tensor contributions ==========
511    if (choice==6) then
512      do ispinor=1,nspinor
513        do ia=1,nincat
514          iashift=3*(ia+ia3-2)
515          do ilmn=1,nlmn
516            do iplex=1,cplex
517              enlk=enlk+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor)
518            end do
519            enlj(1:3)=zero
520            do mu=1,3
521              do iplex=1,cplex
522                enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,6+mu,ilmn,ia,ispinor)
523              end do
524            end do
525            fnlk(iashift+1:iashift+3)=fnlk(iashift+1:iashift+3)+two*enlj(1:3)
526            enlj(1:6)=zero
527            do mu=1,6
528              do iplex=1,cplex
529                enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor)
530              end do
531            end do
532            strnlk(1:6)=strnlk(1:6)+two*enlj(1:6)
533            do mub=1,6
534              mushift=6*(mub-1);nushift=(3*natom+6)*(mub-1)
535              do mua=1,6
536                mu=mushift+mua;nu=nushift+mua
537                do iplex=1,cplex
538                  enlout(nu)=enlout(nu)+two* &
539 &                 (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
540 &                 +dgxdtfac(iplex,mua,ilmn,ia,ispinor)*dgxdt(iplex,mub,ilmn,ia,ispinor))
541                end do
542              end do
543              mushift=36+3*(mub-1);nushift=6+iashift+(3*natom+6)*(mub-1)
544              do mua=1,3
545                mu=mushift+mua;nu=nushift+mua
546                do iplex=1,cplex
547                  enlout(nu)=enlout(nu)+two* &
548 &                 (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
549 &                 +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,6+mua,ilmn,ia,ispinor))
550                end do
551              end do
552            end do
553          end do
554        end do
555      end do
556    end if
557 
558 !  ======== Accumulate the contributions of 2nd-derivatives of E wrt to k ==========
559    if (choice==8) then
560      ABI_ALLOCATE(cft,(3,nlmn))
561      ABI_ALLOCATE(cfu,(3,nlmn))
562      do ispinor=1,nspinor
563        do ia=1,nincat
564          enlj(1:6)=zero
565          do ilmn=1,nlmn
566            do mu=1,6
567              do iplex=1,cplex
568                enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)
569              end do
570            end do
571          end do
572 !        If cplex=1, dgxdt is pure imaginary;
573 !        If cplex_fac=1, dgxdtfac is pure imaginary;
574          if(cplex==2)then
575            cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
576          else
577            cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
578          end if
579          if(cplex_fac==2)then
580            cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor))
581          else
582            cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor))
583          end if
584          do ilmn=1,nlmn
585            do mu=1,6
586              mua=alpha(mu);mub=beta(mu)
587              enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn))
588            end do
589          end do
590          enlout(1:6)=enlout(1:6)+two*enlj(1:6)
591        end do
592      end do
593      ABI_DEALLOCATE(cft)
594      ABI_DEALLOCATE(cfu)
595    end if
596 
597 !  ======== Accumulate the contributions of partial 2nd-derivatives of E wrt to k ==========
598 !  Full derivative wrt to k1, right derivative wrt to k2
599    if (choice==81) then
600      ABI_ALLOCATE(cft,(3,nlmn))
601      ABI_ALLOCATE(cfu,(6,nlmn))
602      ABI_ALLOCATE(enljj,(18))
603      do ispinor=1,nspinor
604        do ia=1,nincat
605          enljj(1:18)=zero
606          if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real
607            cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),gxfac(2,1:nlmn,ia,ispinor))
608          else
609            cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),zero)
610          end if
611          if(cplex==2)then !If cplex=1, d2gxdt is pure real
612            cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor))
613          else
614            cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero)
615          end if
616          do ilmn=1,nlmn
617            do mu=1,3
618              do nu=1,3
619                muu=3*(mu-1)+nu ; mut=gamma(mu,nu)
620                enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn))
621                enljj(2*muu  )=enljj(2*muu  )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn))
622              end do
623            end do
624          end do
625          if(cplex==2)then !If cplex=1, dgxdt is pure imaginary
626            cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
627          else
628            cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
629          end if
630          if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary
631            cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor))
632          else
633            cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor))
634          end if
635          do ilmn=1,nlmn
636            do mu=1,3
637              do nu=1,3
638                muu=3*(mu-1)+nu
639                enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
640                enljj(2*muu  )=enljj(2*muu  )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
641              end do
642            end do
643          end do
644          enlout(1:18)=enlout(1:18)+enljj(1:18)
645        end do
646      end do
647      ABI_DEALLOCATE(cft)
648      ABI_DEALLOCATE(cfu)
649      ABI_DEALLOCATE(enljj)
650    end if
651 
652  end if
653 
654  if (paw_opt==3) then
655 
656 !  ============== Accumulate contribution to <c|S|c> ===============
657    if (choice==1) then
658      do ispinor=1,nspinor
659        do ia=1,nincat
660          do ilmn=1,nlmn
661            do iplex=1,cplex
662              enlout(1)=enlout(1)+gxfac_sij(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor)
663            end do
664          end do
665        end do
666      end do
667    end if
668 
669 !  ============== Accumulate contribution to <c|dS/d_atm_pos|c> ===============
670    if (choice==2.or.choice==23) then
671      ishift=0;if (choice==23) ishift=6
672      do ispinor=1,nspinor
673        do ia=1,nincat
674          enlj(1:3)=zero
675          iashift=3*(ia+ia3-2)
676          do ilmn=1,nlmn
677            do mu=1,3
678              do iplex=1,cplex
679                enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
680              end do
681            end do
682          end do
683          enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3)
684        end do
685      end do
686    end if
687 
688 !  ============== Accumulate contribution to <c|dS/d_strain|c> ===============
689    if (choice==3.or.choice==23) then
690      enlj(1:6)=zero
691      do ispinor=1,nspinor
692        do ia=1,nincat
693          do ilmn=1,nlmn
694            gxfacj(1:cplex)=gxfac_sij(1:cplex,ilmn,ia,ispinor)
695            do iplex=1,cplex
696              enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor)
697            end do
698            do mu=1,6
699              do iplex=1,cplex
700                enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor)
701              end do
702            end do
703          end do
704        end do
705      end do
706      enlout(1:6)=enlout(1:6)+two*enlj(1:6)
707    end if
708 
709 !  ======== Accumulate the contributions of derivatives of <c|S|c> wrt to k ==========
710    if (choice==5) then
711      enlj(1:3)=zero
712 !    If cplex=1, gxfac is real and dgxdt is pure imaginary; thus there is no contribution
713      if(cplex==2)then
714        do ispinor=1,nspinor
715          do ia=1,nincat
716            do ilmn=1,nlmn
717              do mu=1,3
718                do iplex=1,cplex
719                  enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor)
720                end do
721              end do
722            end do
723          end do
724        end do
725      end if
726      enlout(1:3)=enlout(1:3)+two*enlj(1:3)
727    end if
728 
729 !  ====== Accumulate the contributions of left or right derivatives of <c|S|c> wrt to k ==========
730 !  Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k
731    if (choice==51.or.choice==52) then
732      enlj(1:6)=zero
733      do ispinor=1,nspinor
734        do ia=1,nincat
735          if(cplex==2)then
736            do ilmn=1,nlmn
737              do mu=1,3
738                enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) &
739 &               +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor)
740                enlj(2*mu  )=enlj(2*mu  )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) &
741 &               -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
742              end do
743            end do
744          else if (cplex_fac==2) then
745            do ilmn=1,nlmn
746              do mu=1,3
747                enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
748                enlj(2*mu  )=enlj(2*mu  )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
749              end do
750            end do
751          else if (cplex_fac==1) then
752            do ilmn=1,nlmn
753              do mu=1,3
754                enlj(2*mu  )=enlj(2*mu  )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
755              end do
756            end do
757          end if
758        end do
759      end do
760      if (choice==52) then
761        enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6)
762      end if
763      enlout(1:6)=enlout(1:6)+enlj(1:6)
764    end if
765 
766 !  ====== Accumulate contribution to <c|d2S/d_atm_pos d_left_k|c> =========
767    if (choice==54) then
768      ABI_ALLOCATE(enljj,(18))
769      do ispinor=1,nspinor
770        do ia=1,nincat
771          enljj(1:18)=zero
772          iashift=18*(ia+ia3-2)
773          if(cplex==2) then
774            do ilmn=1,nlmn
775              mu=1;nu=1
776              do mua=1,3 ! atm. pos
777                do mub=1,3 ! k
778                  enljj(nu)=enljj(nu) &
779 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) &
780 &                 +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) &
781 &                 +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) &
782 &                 +gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor)
783 
784                  enljj(nu+1)=enljj(nu+1) &
785 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) &
786 &                 -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) &
787 &                 +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) &
788 &                 -gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
789 
790                  mu=mu+1;nu=nu+2
791                end do
792              end do
793            end do
794 !        If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k
795          else
796            do ilmn=1,nlmn
797              mu=1;nu=1
798              do mua=1,3 ! atm. pos
799                do mub=1,3 ! k
800                  enljj(nu+1)=enljj(nu+1) &
801 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) &
802 &                 +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
803                  mu=mu+1;nu=nu+2
804                end do
805              end do
806            end do
807          end if
808          enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18)
809        end do
810      end do
811      ABI_DEALLOCATE(enljj)
812    end if
813 
814 !  ====== Accumulate contribution to <c|d2S/d_dstrain d_right_k|c> =========
815    if (choice==55) then
816      ABI_ALLOCATE(enljj,(36))
817      do ispinor=1,nspinor
818        do ia=1,nincat
819          enljj(1:36)=zero;enlj(:)=zero
820 !        If cplex=1, dgxdt is real for strain, pure imaginary for k;
821 !        If cplex_fac=1, dgxdtfac is pure imaginary for k;
822          if(cplex==2.and.cplex_fac==2) then
823            do ilmn=1,nlmn
824 !            First compute 2nd-derivative contribution
825              mu=1
826              do mua=1,6 ! strain (lambda,nu)
827                mua1=alpha(mua) ! (nu)
828                mua2=beta(mua)  ! (lambda)
829                do mub=1,3 ! k (mu)
830                  muu=3*(gamma(mua1,mub)-1)+mua2
831                  mut=3*(gamma(mua2,mub)-1)+mua1
832                  d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) &
833 &                 +d2gxdt(1:cplex,mut,ilmn,ia,ispinor))
834                  enljj(mu)=enljj(mu) &
835 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) &
836 &                 +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) &
837 &                 +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1)+gxfac_sij(2,ilmn,ia,ispinor)*d2gx(2)
838                  enljj(mu+1)=enljj(mu+1) &
839 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) &
840 &                 -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) &
841 &                 +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(2)-gxfac_sij(2,ilmn,ia,ispinor)*d2gx(1)
842                  mu=mu+2
843                end do
844              end do
845 !            Then store 1st-derivative contribution
846              mu=1
847              do nu=1,3
848                enlj(mu  )=enlj(mu  )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) &
849 &               +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor)
850                enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) &
851 &               -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
852                mu=mu+2
853              end do
854            end do
855 !        If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k
856          else
857            do ilmn=1,nlmn
858              mu=1
859              do mua=1,6 ! strain (lambda,nu)
860                mua1=alpha(mua) ! (nu)
861                mua2=beta(mua)  ! (lambda)
862                do mub=1,3 ! k (mu)
863                  muu=3*(gamma(mua1,mub)-1)+mua2
864                  mut=3*(gamma(mua2,mub)-1)+mua1
865                  d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor))
866                  enljj(mu+1)=enljj(mu+1) &
867 &                 +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) &
868 &                 +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1)
869                  mu=mu+2
870                end do
871              end do
872 !            Then store 1st-derivative contribution
873              mu=1
874              do nu=1,3
875                enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
876                mu=mu+2
877              end do
878            end do
879          end if
880          enlout(1:36)=enlout(1:36)+enljj(1:36)
881          ddkk(1:6)=ddkk(1:6)+enlj(1:6)
882        end do
883      end do
884      ABI_DEALLOCATE(enljj)
885    end if
886 
887 !  ======  Accumulate contribution to <c|d2S/d_k d_k|c> =========
888    if (choice==8) then
889      ABI_ALLOCATE(cft,(3,nlmn))
890      ABI_ALLOCATE(cfu,(3,nlmn))
891      do ispinor=1,nspinor
892        do ia=1,nincat
893          enlj(1:6)=zero
894          do ilmn=1,nlmn
895            do mu=1,6
896              do iplex=1,cplex
897                enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)
898              end do
899            end do
900          end do
901 !        If cplex=1, dgxdt is pure imaginary, dgxdtfac_sij is pure imaginary;
902          if(cplex==2)then
903            cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
904            cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor))
905          else
906            cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
907            cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor))
908          end if
909          do ilmn=1,nlmn
910            do mu=1,6
911              mua=alpha(mu);mub=beta(mu)
912              enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn))
913            end do
914          end do
915          enlout(1:6)=enlout(1:6)+two*enlj(1:6)
916        end do
917      end do
918      ABI_DEALLOCATE(cft)
919      ABI_DEALLOCATE(cfu)
920    end if
921 
922 !  ======  Accumulate contribution to <c|d/d_k[d(right)S/d_k]|c> =========
923 !  Full derivative wrt to k1, right derivative wrt to k2
924    if (choice==81) then
925      ABI_ALLOCATE(cft,(3,nlmn))
926      ABI_ALLOCATE(cfu,(6,nlmn))
927      ABI_ALLOCATE(enljj,(18))
928      do ispinor=1,nspinor
929        do ia=1,nincat
930          enljj(1:18)=zero
931          if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real
932            cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),gxfac_sij(2,1:nlmn,ia,ispinor))
933          else
934            cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),zero)
935          end if
936          if(cplex==2)then !If cplex=1, d2gxdt is pure real
937            cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor))
938          else
939            cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero)
940          end if
941          do ilmn=1,nlmn
942            do mu=1,3
943              do nu=1,3
944                muu=3*(mu-1)+nu ; mut=gamma(mu,nu)
945                enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn))
946                enljj(2*muu  )=enljj(2*muu  )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn))
947              end do
948            end do
949          end do
950          if(cplex==2)then !If cplex=1, dgxdt is pure imaginary
951            cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
952          else
953            cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
954          end if
955          if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary
956            cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor))
957          else
958            cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor))
959          end if
960          do ilmn=1,nlmn
961            do mu=1,3
962              do nu=1,3
963                muu=3*(mu-1)+nu
964                enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
965                enljj(2*muu  )=enljj(2*muu  )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
966              end do
967            end do
968          end do
969          enlout(1:18)=enlout(1:18)+enljj(1:18)
970        end do
971      end do
972      ABI_DEALLOCATE(cft)
973      ABI_DEALLOCATE(cfu)
974      ABI_DEALLOCATE(enljj)
975    end if
976 
977  end if
978 
979 end subroutine opernld_ylm