TABLE OF CONTENTS


ABINIT/ddkten [ Functions ]

[ Top ] [ Functions ]

NAME

 ddkten

FUNCTION

 Compact or decompact the tensors related to the ffnl(:,1,...)
 part of the ddk operator, taking into account the direction
 of the ddk perturbation.

INPUTS

  compact= if 1, compact from tmpfac
  idir=direction of the ddk perturbation
  rank=0,1,2, or 3 = rank of tmpfac tensor, also angular momentum (=l)

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output:
  temp(2,(rank*(rank+1))/2)=compacted tensor
    for l=1, just a scalar
    for l=2, a vector
  tmpfac(2,(rank+1)*(rank+2)/2)=decompacted tensor
    for l=1, a vector
    for l=2, a symmetric matrix, stored as
     (1 . .)
     (6 2 .)
     (5 4 3)

NOTES

 For l=0, there is no contribution.

SOURCE

1536 subroutine ddkten(compact,idir,rank,temp,tmpfac)
1537 
1538 !Arguments ------------------------------------
1539 !scalars
1540  integer,intent(in) :: compact,idir,rank
1541 !arrays
1542  real(dp),intent(inout) :: temp(2,(rank*(rank+1))/2)
1543  real(dp),intent(inout) :: tmpfac(2,((rank+1)*(rank+2))/2)
1544 
1545 !Local variables-------------------------------
1546 !scalars
1547  character(len=500) :: msg
1548 
1549 ! *************************************************************************
1550 
1551  if(rank/=1 .and. rank/=2 .and. rank/=3)then
1552    write(msg, '(a,i10,a,a,a)' )&
1553    'Input rank=',rank,' not allowed.',ch10,&
1554    'Possible values are 1,2,3 only.'
1555    ABI_BUG(msg)
1556  end if
1557 
1558 !Take care of p angular momentum
1559  if(rank==1)then
1560 
1561 !  Compaction tmpfac -> temp
1562    if(compact==1)then
1563      temp(:,1)=tmpfac(:,idir)
1564 
1565 !    Decompaction temp -> tmpfac
1566    else
1567      tmpfac(:,1:3)=0.0d0
1568      tmpfac(:,idir)=temp(:,1)
1569    end if
1570 
1571 !  Take care of d angular momentum
1572 !  rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1573 
1574  else if(rank==2)then
1575 
1576 !  Compaction tmpfac -> temp
1577    if(compact==1)then
1578      if(idir==1)then
1579 !      Count the number of non-zero derivatives with respect to k(idir)
1580 !      The factor of 2 on the diagonal comes from the derivative with
1581 !      respect to the first K then to the second K
1582        temp(:,1)=2.0d0*tmpfac(:,1); temp(:,2)=tmpfac(:,6); temp(:,3)=tmpfac(:,5)
1583      else if(idir==2)then
1584        temp(:,2)=2.0d0*tmpfac(:,2); temp(:,1)=tmpfac(:,6); temp(:,3)=tmpfac(:,4)
1585      else if(idir==3)then
1586        temp(:,3)=2.0d0*tmpfac(:,3); temp(:,1)=tmpfac(:,5); temp(:,2)=tmpfac(:,4)
1587      end if
1588 !    Decompaction temp -> tmpfac
1589    else
1590      tmpfac(:,1:6)=0.0d0
1591      tmpfac(:,idir)=2.0d0*temp(:,idir)
1592      if(idir==1)then
1593        tmpfac(:,5)=temp(:,3); tmpfac(:,6)=temp(:,2)
1594      else if(idir==2)then
1595        tmpfac(:,4)=temp(:,3); tmpfac(:,6)=temp(:,1)
1596      else if(idir==3)then
1597        tmpfac(:,4)=temp(:,2); tmpfac(:,5)=temp(:,1)
1598      end if
1599    end if
1600 
1601 !  Take care of f angular momentum
1602  else if(rank==3)then
1603 !  rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
1604 !  rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1605 
1606 !  Compaction tmpfac -> temp
1607    if(compact==1)then
1608      if(idir==1)then
1609 !      Count the number of non-zero derivatives with respect to k(idir)
1610        temp(:,1)=3.0d0*tmpfac(:,1)
1611        temp(:,2:4)=tmpfac(:,2:4)
1612        temp(:,5:6)=2.0d0*tmpfac(:,5:6)
1613      else if(idir==2)then
1614        temp(:,6)=2.0d0*tmpfac(:,2)
1615        temp(:,4)=2.0d0*tmpfac(:,9)
1616        temp(:,5)=tmpfac(:,4)
1617        temp(:,1)=tmpfac(:,6)
1618        temp(:,3)=tmpfac(:,8)
1619        temp(:,2)=3.0d0*tmpfac(:,7)
1620      else if(idir==3)then
1621        temp(:,3)=3.0d0*tmpfac(:,10)
1622        temp(:,5)=2.0d0*tmpfac(:,3)
1623        temp(:,4)=2.0d0*tmpfac(:,8)
1624        temp(:,6)=tmpfac(:,4)
1625        temp(:,1)=tmpfac(:,5)
1626        temp(:,2)=tmpfac(:,9)
1627      end if
1628 !    Decompaction temp -> tmpfac
1629    else
1630      tmpfac(:,1:10)=0.0d0
1631      if(idir==1)then
1632        tmpfac(:,1)=3.0d0*temp(:,1)
1633        tmpfac(:,2:4)=temp(:,2:4)
1634        tmpfac(:,5:6)=2.0d0*temp(:,5:6)
1635      else if(idir==2)then
1636        tmpfac(:,2)=2.0d0*temp(:,6)
1637        tmpfac(:,9)=2.0d0*temp(:,4)
1638        tmpfac(:,4)=temp(:,5)
1639        tmpfac(:,6)=temp(:,1)
1640        tmpfac(:,8)=temp(:,3)
1641        tmpfac(:,7)=3.0d0*temp(:,2)
1642      else if(idir==3)then
1643        tmpfac(:,10)=3.0d0*temp(:,3)
1644        tmpfac(:,3)=2.0d0*temp(:,5)
1645        tmpfac(:,8)=2.0d0*temp(:,4)
1646        tmpfac(:,4)=temp(:,6)
1647        tmpfac(:,5)=temp(:,1)
1648        tmpfac(:,9)=temp(:,2)
1649      end if
1650    end if
1651 
1652  end if
1653 
1654 end subroutine ddkten

ABINIT/m_nonlop_pl [ Modules ]

[ Top ] [ Modules ]

NAME

  nonlop_pl

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, GZ, MT, FF, 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

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module m_nonlop_pl
22 
23  use defs_basis
24  use m_errors
25  use m_abicore
26  use m_xmpi
27  use m_contistr01
28  use m_contistr03
29  use m_contistr12
30  use m_contstr21
31  use m_contstr23
32  use m_contstr25
33  use m_contstr25a
34  use m_contstr26
35  use m_metstr
36  use m_opernl
37 
38  use defs_abitypes, only : MPI_type
39  use m_geometry,   only : strconv
40  use m_kg,         only : ph1d3d
41  use m_contract,   only : cont22cso, cont22so, cont24, cont33cso, cont33so, cont35, cont22, cont3, cont13, &
42                           metcon, metcon_so, metric_so
43  implicit none
44 
45  private

ABINIT/nonlop_pl [ Functions ]

[ Top ] [ Functions ]

NAME

 nonlop_pl

FUNCTION

 * Compute application of a nonlocal operator Vnl in order to get:
    - contracted elements (energy, forces, stresses, ...), if signs=1
    - a function in reciprocal space (|out> = Vnl|in>), if signs=2
   Operator Vnl, as the following form:
    $Vnl=sum_{R,lmn,l''m''n''} {|P_{Rlmn}> Enl^{R}_{lmn,l''m''n''} <P_{Rl''m''n''}|}$
   Operator Vnl is -- in the typical case -- the nonlocal potential.
   - With norm-conserving pseudopots, $Enl^{R}_{lmn,l''m''n''}$ is the
     Kleinmann-Bylander energy $Ekb^{R}_{ln}$.
   - The |P_{Rlmn}> are the projector functions.
 * This routine uses Legendre polynomials Pl to express Vnl.

INPUTS

  choice: chooses possible output:
    choice=1 => a non-local energy contribution
          =2 => a gradient with respect to atomic position(s)
          =3 => a gradient with respect to strain(s)
          =23=> a gradient with respect to atm. pos. and strain(s)
          =4 => a gradient and 2nd derivative with respect to atomic pos.
          =5 => a gradient with respect to k wavevector
          =6 => 2nd derivatives with respect to strain
  dimekb1,dimekb2=dimensions of ekb (see ekb)
  dimffnlin=second dimension of ffnlin (1+number of derivatives)
  dimffnlout=second dimension of ffnlout (1+number of derivatives)
  ekb(dimekb1,dimekb2,nspinortot**2)= (Real) Kleinman-Bylander energies (hartree)
                                   dimekb1=lmnmax  -  dimekp2=ntypat
  ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used
          for the application of the nonlocal operator to the |in> vector
  ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used
          for the application of the nonlocal operator to the |out> vector
  gmet(3,3)=metric tensor for G vecs (in bohr**-2)
  gprimd(3,3)=dimensional reciprocal space primitive translations
   (bohr^-1)
  idir=direction of the - atom to be moved in the case (choice=2,signs=2),
                        - k point direction in the case (choice=5,signs=2)
                        - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
  indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
  istwf_k=option parameter that describes the storage of wfs
  kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector
  kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector
  kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector
  kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector
  kptin(3)=k point in terms of recip. translations, for the |in> vector
  kptout(3)=k point in terms of recip. translations, for the |out> vector
  lmnmax=max. number of (l,m,n) components over all types of atoms
  matblk=dimension of the arrays ph3din and ph3dout
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nattyp(ntypat)=number of atoms of each type
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpgin,nkpgout=second sizes of arrays kpgin/kpgout
  nloalg(3)=governs the choice of the algorithm for nonlocal operator
  nnlout=dimension of enlout: choice=1=>nnlout=1   choice=2=>nnlout=3*natom
                              choice=3=>nnlout=6   choice=4=>nnlout=6*natom
                              choice=5=>nnlout=1   choice=6=>nnlout=6*(3*natom+6)
                              choice=23=>nnlout=6+3*natom
  npwin=number of planewaves for given k point, for the |in> vector
  npwout=number of planewaves for given k point, for the |out> vector
  nspinor=number of spinorial components of the wavefunctions on current proc
  nspinortot=total number of spinorial components of the wavefunctions
  ntypat=number of types of atoms in cell
  only_SO=flag to calculate only the SO part in nonlop
  phkxredin(2,natom)=phase factors exp(2 pi kptin.xred)
  phkxredout(2,natom)=phase factors exp(2 pi kptout.xred)
  ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
  ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in)
  ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out)
  --- pspso removed in beautification because it was unused ---
  pspso(ntypat)=spin-orbit characteristic for each atom type
  -------------------------------------------------------------
  signs= if 1, get contracted elements (energy, forces, stress, ...)
         if 2, applies the non-local operator to a function in reciprocal space
  ucvol=unit cell volume (bohr^3)
  vectin(2,nspinor*npwin)=input cmplx wavefunction coefficients <G|Cnk>

OUTPUT

  ==== if (signs==1) ====
     enlout(nnlout)= contribution of this state to the nl part
                     of the following properties:
       if choice=1 : enlout(1)               -> the energy
       if choice=2 : enlout(1:3*natom)       -> the forces
       if choice=3 : enlout(1:6)             -> the stresses
       if choice=23: enlout(1:6+3*natom)     -> the forces and the stresses
       if choice=4 : enlout(1:6*natom)       -> the frozen wf part of dynam. matrix
       if choice=6 : enlout(1:6*(3*natom+6)) -> the frozen wf part of elastic tensor
  ==== if (signs==2) ====
     vectout(2,nspinor*npwout)= result of the aplication of the nl operator
                                or one of its derivative to the input vect.:
       if choice=1 : Vnl |vectin>
       if choice=2 : dVnl/d(xred(idir,iatom) |vectin> (xred=reduced atm. pos.)
       if choice=3 : dVnl/d(strain(idir)) |vectin>    (symmetric strain =>idir=1...6)
       if choice=5 : dVnl/dk(idir) |vectin>           (k wavevector)

NOTES

 In the case signs=1, the array vectout is not used, nor modified
 so that the same array as vectin can be used as a dummy argument;
 the same is true for the pairs npwin-npwout, ffnlin-ffnlout,
 kgin-kgout, ph3din-ph3dout, phkredin-phkxredout).

 Calculation includes contributions to grads of Etot wrt coord and
 wrt strains for l=0,1,2,3.

WARNINGS

  - Warning 1: This routine is in a transient state, during the
    time of the implementation of the spin-orbit coupling...
    In particular, the OMP parallelisation is still missing,
    but it matters here only when nspinor==2.
  - Warning 2: the order of atoms is governed by atindx

SOURCE

 170 subroutine nonlop_pl(choice,dimekb1,dimekb2,dimffnlin,dimffnlout,ekb,enlout,&
 171 &                     ffnlin,ffnlout,gmet,gprimd,idir,indlmn,istwf_k,kgin,kgout,kpgin,kpgout,&
 172 &                     kptin,kptout,lmnmax,matblk,mgfft,mpi_enreg,mpsang,mpssoang,&
 173 &                     natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,npwin,npwout,nspinor,nspinortot,&
 174 &                     ntypat,only_SO,phkxredin,phkxredout,ph1d,ph3din,ph3dout,signs,&
 175 &                     ucvol,vectin,vectout)
 176 
 177 !Arguments ------------------------------------
 178 !This type is defined in defs_mpi
 179 !The (inout) classification below is misleading; mpi_enreg is temporarily
 180 ! changed but reset to its initial condition before exiting.
 181 !scalars
 182  integer,intent(in) :: choice,dimekb1,dimekb2,dimffnlin,dimffnlout,idir,istwf_k
 183  integer,intent(in) :: lmnmax,matblk,mgfft,mpsang,mpssoang,natom,nkpgin,nkpgout
 184  integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,only_SO,signs
 185  real(dp),intent(in) :: ucvol
 186  type(MPI_type),intent(in) :: mpi_enreg
 187 !arrays
 188  integer,intent(in) :: indlmn(6,lmnmax,ntypat),kgin(3,npwin),kgout(3,npwout)
 189  integer,intent(in) :: nattyp(ntypat),ngfft(18),nloalg(3) !,pspso(ntypat) UNUSED
 190  real(dp),intent(in) :: ekb(dimekb1,dimekb2,nspinortot**2)
 191  real(dp),intent(in) :: ffnlin(npwin,dimffnlin,lmnmax,ntypat)
 192  real(dp),intent(in) :: ffnlout(npwout,dimffnlout,lmnmax,ntypat),gmet(3,3)
 193  real(dp),intent(in) :: gprimd(3,3),kpgin(npwin,nkpgin),kpgout(npwout,nkpgout)
 194 !real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom) !vz_d
 195  real(dp),intent(in) :: kptin(3),kptout(3) !vz_d
 196  real(dp),intent(in) :: ph1d(2,*) !vz_d
 197  real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom)
 198  real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk)
 199  real(dp),intent(inout) :: vectin(:,:)
 200  real(dp),intent(out) :: enlout(:) !vz_i
 201  real(dp),intent(inout) :: vectout(:,:) !vz_i
 202 
 203 !Local variables-------------------------------
 204 !mlang is the maximum number of different angular momenta
 205 !(mlang=4 means s,p,d,f)
 206 ! Note : in a future version, one should adjust mlang to mpsang.
 207 !mlang2 is the maximum number of unique tensor components for a tensor
 208 !of rank (mlang-1) with index range 1-3
 209 !mlang3 is the maximum number of unique tensor components summed over
 210 !all tensors of rank 0 through mlang-1.
 211 !mlang4 is the total number of additional unique tensor components
 212 !related to strain gradients, ranks 2 through mlang+1.
 213 !mlang6 is the total number of additional unique tensor components
 214 !related to strain 2nd derivaives, ranks 4 through mlang+3.
 215 !mlang1 is the total number of certain additional unique tensor components
 216 !related to internal strain, ranks 1 through mlang
 217 !mlang5 is the total number of other additional unique tensor components
 218 !related to internal strain, ranks 1 through mlang
 219 !scalars
 220  integer,parameter :: mlang=4
 221 ! MG: I tried to use parameters instead of saved variables but [tutorespfn][trf2_1] gets stuck on milou_g95_snofbfpe
 222  integer,save :: mlang1=((mlang+1)*(mlang+2)*(mlang+3))/6-1
 223  !integer,save :: mlang2=(mlang*(mlang+1))/2 ! Unused
 224  integer,save :: mlang3=(mlang*(mlang+1)*(mlang+2))/6
 225  integer,save :: mlang4=((mlang+2)*(mlang+3)*(mlang+4))/6-4
 226  integer,save :: mlang5=((mlang+3)*(mlang+4)*(mlang+5))/6-10
 227  integer,save :: mlang6=((mlang+4)*(mlang+5)*(mlang+6))/6-20
 228  integer :: compact,ia,ia1,ia2,ia3,ia4,ia5,ierr,iest,ig,ii,ilang,ilang2,ilmn
 229  integer :: iln,iln0,indx,iproj,ipsang,ishift,isp,ispin,ispinor,ispinor_index,ispinp
 230  integer :: istr,istr1,istr2,iterm,itypat,jj,jjk,jjs,jjs1,jjs2,jjs3,jjs4,jjstr,jspin
 231  integer :: mincat,mproj,mu,mumax,n1,n2,n3,ndgxdt,ndgxdtfac,nincat,nlang
 232  integer :: nproj,nspinso,rank
 233  integer :: sign,spaceComm,  isft
 234  real(dp) :: e2nl,e2nldd,enlk
 235  character(len=500) :: msg
 236 !arrays
 237  integer,allocatable :: indlmn_s(:,:,:),jproj(:)
 238  real(dp) :: amet(2,3,3,2,2),amet_lo(3,3),e2nl_tmp(6),eisnl(3),rank2(6)
 239  real(dp) :: rank2c(2,6),strsnl(6),strsnl_out(6),strsso(6,3),strssoc(6),trace(2)!,tsec(2)
 240  real(dp),allocatable :: d2gxdis(:,:,:,:,:),d2gxdis_s(:,:,:,:)
 241  real(dp),allocatable :: d2gxds2(:,:,:,:,:),d2gxds2_s(:,:,:,:)
 242  real(dp),allocatable :: dgxdis(:,:,:,:,:),dgxdis_s(:,:,:,:),dgxds(:,:,:,:,:)
 243  real(dp),allocatable :: dgxds_s(:,:,:,:),dgxdsfac(:,:,:,:,:)
 244  real(dp),allocatable :: dgxdt(:,:,:,:,:,:),dgxdt_s(:,:,:,:,:)
 245  real(dp),allocatable :: dgxdtfac(:,:,:,:,:),ekb_s(:,:),gxa(:,:,:,:,:)
 246  real(dp),allocatable :: gxa_s(:,:,:,:),gxafac(:,:,:,:),pauli(:,:,:,:)
 247  real(dp),allocatable :: temp(:,:),tmpfac(:,:),vectin_s(:,:),vectout_s(:,:)
 248  real(dp),allocatable :: wt(:,:)
 249 
 250 ! **********************************************************************
 251 
 252  ABI_UNUSED(mgfft)
 253 
 254 !Test: spin orbit not allowed for choice=5,6
 255  if (nspinortot==2 .and. choice==6 ) then
 256    ABI_BUG('For nspinortot=2, choice=6 is not yet allowed.')
 257  end if
 258 
 259  if ((choice<1 .or. choice>6) .and. choice/=23 ) then
 260    write(msg,'(a,i0)')'Does not presently support this choice=',choice
 261    ABI_BUG(msg)
 262  end if
 263 
 264 !Test: choice 51 and 52 only allowed with nonlop_ylm
 265 !JWZ, 01-Sep-08
 266  if (choice==51 .or. choice==52) then
 267    ABI_BUG('choice 51 or 52 is not yet allowed.')
 268  end if
 269 
 270 !Define dimension of work arrays.
 271  mincat=min(NLO_MINCAT,maxval(nattyp))
 272  mproj=maxval(indlmn(3,:,:))
 273  ABI_MALLOC(temp,(2,mlang4))
 274  ABI_MALLOC(tmpfac,(2,mlang4))
 275  ABI_MALLOC(wt,(mlang,mproj))
 276  ABI_MALLOC(jproj,(mlang))
 277  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
 278  ABI_MALLOC(ekb_s,(mlang,mproj))
 279  ABI_MALLOC(indlmn_s,(6,lmnmax,ntypat))
 280 
 281 !Eventually compute the spin-orbit metric tensor:
 282  if (mpssoang>mpsang) then
 283    ABI_MALLOC(pauli,(2,2,2,3))
 284    call metric_so(amet,gprimd,pauli)
 285  end if
 286 
 287 !Allocate array gxa (contains projected scalars).
 288  ABI_MALLOC(gxa,(2,mlang3,mincat,mproj,nspinortot))
 289  if(nspinor==2)  then
 290    ABI_MALLOC(gxa_s,(2,mlang3,mincat,mproj))
 291  else
 292    ABI_MALLOC(gxa_s,(0,0,0,0))
 293  end if
 294 
 295  ABI_MALLOC(gxafac,(2,mlang3,mincat,mproj))
 296  gxa(:,:,:,:,:)=zero
 297 
 298 !If choice==2 : first-order atomic displacements
 299 !If signs==2, only one direction is considered
 300 !If signs==1, the three directions are considered
 301 !If choice==4 and signs==1 : second-order atomic displacements,
 302 !the nine components are considered
 303 !If choice==5 and signs==2 : ddk
 304 !component 1 -> from ffnl(:,2,...)
 305 !component 2 -> from ffnl(:,1,...) (too much space is booked for this
 306 !case, since the number of angular momenta is smaller than mlang3, but it is easier)
 307  ndgxdt=0
 308  if(signs==2 .and. choice==2) ndgxdt=1
 309  if(signs==1 .and. (choice==2.or.choice==23)) ndgxdt=3
 310  if(choice==4) ndgxdt=9
 311  if(choice==5) ndgxdt=2
 312 !Allocate dgxdt (contains derivatives of gxa with respect to atomic displacements or ddk).
 313  ABI_MALLOC(dgxdt,(2,ndgxdt,mlang3,mincat,mproj,nspinortot))
 314  dgxdt(:,:,:,:,:,:)=zero
 315  if(nspinor==2)then
 316    ABI_MALLOC(dgxdt_s,(2,ndgxdt,mlang3,mincat,mproj))
 317    dgxdt_s(:,:,:,:,:)=zero
 318  else
 319    ABI_MALLOC(dgxdt_s,(0,0,0,0,0))
 320  end if
 321  ndgxdtfac=0
 322  if(signs==2 .and. choice==2) ndgxdtfac=1
 323  if(choice==4) ndgxdtfac=3
 324  if(choice==5) ndgxdtfac=2
 325  ABI_MALLOC(dgxdtfac,(2,ndgxdtfac,mlang3,mincat,mproj))
 326 
 327 !Allocate dgxds (contains derivatives of gxa with respect to strains).
 328  ABI_MALLOC(dgxds,(2,mlang4,mincat,mproj,nspinor))
 329  dgxds(:,:,:,:,:)=zero
 330  ABI_MALLOC(dgxdsfac,(2,mlang4,mincat,mproj,nspinor))
 331  if(choice==6) then
 332    ABI_MALLOC(dgxdis,(2,mlang1,mincat,mproj,nspinor))
 333    ABI_MALLOC(d2gxdis,(2,mlang5,mincat,mproj,nspinor))
 334    ABI_MALLOC(d2gxds2,(2,mlang6,mincat,mproj,nspinor))
 335  else
 336    ABI_MALLOC(dgxdis ,(0,0,0,0,0))
 337    ABI_MALLOC(d2gxdis,(0,0,0,0,0))
 338    ABI_MALLOC(d2gxds2,(0,0,0,0,0))
 339  end if
 340  ABI_MALLOC(dgxds_s  ,(0,0,0,0))
 341  ABI_MALLOC(dgxdis_s ,(0,0,0,0))
 342  ABI_MALLOC(d2gxdis_s,(0,0,0,0))
 343  ABI_MALLOC(d2gxds2_s,(0,0,0,0))
 344  if(nspinor==2)then
 345    ABI_FREE(dgxds_s)
 346    ABI_MALLOC(dgxds_s,(2,mlang4,mincat,mproj))
 347    dgxds_s(:,:,:,:)=zero
 348    if(choice==6) then
 349      ABI_FREE(dgxdis_s)
 350      ABI_FREE(d2gxdis_s)
 351      ABI_FREE(d2gxds2_s)
 352      ABI_MALLOC(dgxdis_s,(2,mlang1,mincat,mproj))
 353      ABI_MALLOC(d2gxdis_s,(2,mlang5,mincat,mproj))
 354      ABI_MALLOC(d2gxds2_s,(2,mlang6,mincat,mproj))
 355    else
 356    end if
 357  end if
 358 
 359 !Zero out some arrays
 360  if(choice==2 .or. choice==4 .or. choice==5 .or. choice==6 .or. choice==23) enlout(:)=0.0d0
 361 
 362  if(signs==2) vectout(:,:)=zero
 363 
 364  !if(choice==3.or.choice==23) then
 365  !  strsnl(:)=zero
 366  !  if(mpssoang>mpsang) strsso(:,:)=zero
 367  !end if
 368  enlk=zero
 369  strsso = zero
 370  strsnl = zero
 371 
 372 !In the case vectin is a spinor, split its second part.
 373 !Also, eventually take into account the storage format of the wavefunction
 374 !(the original vector will be restored before leaving the routine,
 375 !except for the vectin(2,1) component with istwf_k==2,
 376 !that should vanish)
 377 !In sequential, treat the second spinor part first
 378  if (nspinor==2)then
 379    ABI_MALLOC(vectin_s,(2,npwin))
 380    ABI_MALLOC(vectout_s,(2,npwout))
 381 
 382    isft = npwin;if (mpi_enreg%nproc_spinor>1) isft=0
 383 
 384 !  Initialize it
 385 !$OMP PARALLEL DO
 386    do ig=1,npwin
 387      vectin_s(1,ig)=vectin(1,ig+isft)
 388      vectin_s(2,ig)=vectin(2,ig+isft)
 389    end do
 390 
 391 !  Take into account the storage
 392    if(istwf_k/=1)then
 393      call scalewf_nonlop(istwf_k,mpi_enreg,npwin,1,vectin_s)
 394    end if
 395  end if ! nspinortot==2
 396 
 397 !Treat the first spinor part now
 398  if(istwf_k/=1) then
 399    call scalewf_nonlop(istwf_k,mpi_enreg,npwin,1,vectin)
 400  end if
 401 
 402 
 403 !Big loop on atom types.
 404  ia1=1
 405  do itypat=1,ntypat
 406 
 407 !  Get atom loop indices for different types:
 408    ia2=ia1+nattyp(itypat)-1
 409 
 410 !  Cut the sum on different atoms in blocks, to allow memory saving.
 411 !  Inner summations on atoms will be done from ia3 to ia4.
 412 !  Note that the maximum range from ia3 to ia4 is mincat (maximum
 413 !  increment of atoms).
 414    do ia3=ia1,ia2,mincat
 415      ia4=min(ia2,ia3+mincat-1)
 416 !    Give the increment of number of atoms in this subset.
 417      nincat=ia4-ia3+1
 418 
 419 !    Prepare the phase factors for the atoms between ia3 and ia4 :
 420 !    For nloalg(2)<=0, they were not prepared previously,it is needed to
 421 !    compute them again.
 422      if(nloalg(2)<=0)then
 423 !      For nloalg(2)==0, it is needed to compute the phase factors.
 424        if(mincat>matblk)then
 425          write(msg,'(a,a,a,i4,a,i4,a)')&
 426           'With nloc_mem<=0, mincat must be less than matblk.',ch10,&
 427           'Their value is ',mincat,' and ',matblk,'.'
 428          ABI_BUG(msg)
 429        end if
 430        call ph1d3d(ia3,ia4,kgin,matblk,natom,npwin,n1,n2,n3,phkxredin,ph1d,ph3din)
 431      end if
 432 
 433 !    Here begins the different treatment for the scalar-relativistic
 434 !    part(s) and the spin-orbit part.
 435 !    Loop on ispinor : 1 for scalar-Relativistic, 2 for spin-orbit
 436      nspinso=1;if (mpssoang>mpsang) nspinso=2
 437 
 438      ! Change nspinso if collinear run or if nspinor == 2 and SOC is not wanted.
 439      ! TODO: The last check requires pspso
 440      if (nspinortot == 1) nspinso = 1
 441 
 442      do ispinor=1,nspinso
 443        ispinor_index=ispinor
 444        if (mpi_enreg%paral_spinor==1) ispinor_index=mpi_enreg%me_spinor+1
 445 
 446 !
 447 !      mjv 25 6 2008: if only_SO == 1 or 2 skip the scalar relativistic terms
 448 !      only output the spin orbit ones
 449 !
 450        if (ispinor==1 .and. only_SO>0) cycle
 451 
 452 !      Select scalar-relativistic or spin-orbit KB-energies:
 453        ekb_s(:,:)=zero ; wt(:,:)=zero
 454 !      Loop over (l,n) values (loop over l,m,n and test on l,n)
 455        iln0=0 ; jproj(:)=0 ; nlang=0
 456        indlmn_s(:,:,itypat)=0
 457        do ilmn=1,lmnmax
 458          if(ispinor/=indlmn(6,ilmn,itypat))cycle
 459          indlmn_s(:,ilmn,itypat)=indlmn(:,ilmn,itypat)
 460          iln =indlmn(5,ilmn,itypat)
 461          if (iln>iln0) then
 462            iln0=iln
 463            ipsang=indlmn(1,ilmn,itypat)+1
 464 !          DEBUG
 465 !          write(std_out,*)' nonlop : ipsang,ilmn,itypat,ispinor=',ipsang,ilmn,itypat,ispinor
 466 !          ENDDEBUG
 467            iproj=indlmn(3,ilmn,itypat)
 468 !          This shift is not needed anymore
 469 !          if (ispinor==2) ipsang=indlmn(1,ilmn,itypat)-mpsang+2
 470            ekb_s(ipsang,iproj)=ekb(iln,itypat,ispinor)
 471            wt(ipsang,iproj)=4.d0*pi/ucvol*dble(2*ipsang-1)*ekb_s(ipsang,iproj)
 472 !
 473 !          mjv 27 6 2008: if only_SO == 2 remove the factor of l in the operator
 474 !
 475            if (only_SO == 2) then
 476              wt(ipsang,iproj)=4.d0*pi/ucvol*ekb_s(ipsang,iproj)
 477            end if
 478            jproj(ipsang)=max(jproj(ipsang),iproj)
 479            if(iproj>0)nlang=max(nlang,ipsang)
 480          end if
 481        end do ! ilmn
 482 
 483 
 484 !      If nlang is not 0, then some non-local part is to be applied for that type of atom.
 485        if (nlang/=0) then
 486 !        Operate with the non-local potential on the wavefunction, in order
 487 !        to get projected quantities. Call different routines opernl2,
 488 !        opernl3, opernl4x which corresponds to different writings of the
 489 !        same numerical operations. There is still optimisation left for
 490 !        the case istwf_k/=1 (up to a factor 2 in speed).
 491 !        call timab(74+choice,1,tsec)
 492          sign=1
 493          gxa(:,:,:,:,:)=zero
 494          dgxdt(:,:,:,:,:,:)=zero
 495          dgxds(:,:,:,:,:)=zero
 496 
 497 !        Only the first spinorial component of vectin is taken into account first
 498          ispin=1;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index
 499          if(nloalg(1)==2) then
 500            call opernl2(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
 501 &           ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
 502 &           jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 503 &           mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
 504 &           ntypat,ph3din,sign,vectin)
 505          else if(nloalg(1)==3) then
 506            call opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
 507 &           ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
 508 &           jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 509 &           mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
 510 &           ntypat,ph3din,sign,vectin)
 511          else if(nloalg(1)==4) then
 512            call opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
 513 &           ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
 514 &           jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 515 &           mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
 516 &           ntypat,ph3din,vectin)
 517          end if
 518 !        This duplication of the opernl calls is needed to avoid copying
 519 !        vectin, with a detrimental effect on speed.
 520          if (nspinor==2)then
 521            ispin=2
 522            if(nloalg(1)==2) then
 523              call opernl2(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,&
 524 &             ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
 525 &             jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 526 &             mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
 527 &             ntypat,ph3din,sign,vectin_s)
 528            else if(nloalg(1)==3) then
 529              call opernl3(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,&
 530 &             ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
 531 &             jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 532 &             mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
 533 &             ntypat,ph3din,sign,vectin_s)
 534            else if(nloalg(1)==4) then
 535              call opernl4a(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,&
 536 &             ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
 537 &             jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
 538 &             mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,&
 539 &             ntypat,ph3din,vectin_s)
 540            end if
 541            dgxds(:,:,:,:,ispin)=dgxds_s(:,:,:,:)
 542            dgxdt(:,:,:,:,:,ispin)=dgxdt_s(:,:,:,:,:)
 543            gxa(:,:,:,:,ispin)=gxa_s(:,:,:,:)
 544          end if
 545 
 546 !        Parallelism stuff
 547          spaceComm=mpi_enreg%comm_fft
 548          call xmpi_sum(dgxds,spaceComm,ierr)
 549          if (mpi_enreg%paral_spinor==1) then
 550            spaceComm=mpi_enreg%comm_spinorfft
 551            jspin=3-ispin
 552            gxa(:,:,:,:,ispin)=gxa(:,:,:,:,1)
 553            gxa(:,:,:,:,jspin)=zero
 554            if ( ndgxdt>0)then
 555              dgxdt(:,:,:,:,:,ispin)=dgxdt(:,:,:,:,:,1)
 556              dgxdt(:,:,:,:,:,jspin)=zero
 557            end if
 558          end if
 559 
 560          call xmpi_sum(gxa,spaceComm,ierr)
 561          if(ndgxdt>0) then
 562            call xmpi_sum(dgxdt,spaceComm,ierr)
 563          end if
 564 
 565 !        XG030513 : MPIWF, at this place, one should perform the reduction
 566 !        and spread of data of gxa, dgxdt and dgxds
 567 
 568 !        Loop over spins:
 569          do isp=1,nspinor
 570            ispin=isp;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index
 571 
 572 !          Perform contractions for the various tensors (d)gx?, producing the
 573 !          contracted tensors (d)gx?fac to be passed back to opernl:
 574            do ia=1,nincat
 575              do ilang=1,nlang
 576                nproj=jproj(ilang)
 577                if(nproj/=0) then
 578                  ilang2=(ilang*(ilang+1))/2
 579                  do iproj=1,nproj
 580 
 581 !                  The rank of the tensor gxa equals l:
 582                    rank=ilang-1
 583 !                  jjs gives the starting address of the relevant components
 584                    jjs=1+((ilang-1)*ilang*(ilang+1))/6
 585                    if (ilang>4) then
 586                      write(msg,'(a,i0)')' ilang must fall in range [1..4] but value is ',ilang
 587                      ABI_BUG(msg)
 588                    end if
 589 
 590 !                  Metric & spinorial contraction from gxa to gxafac. The treatment
 591 !                  is different for the scalar-relativistic and spin-orbit parts.
 592                    if(ispinor==1) then
 593 !                    ------ Scalar-Relativistic ------
 594                      temp(:,1:((rank+1)*(rank+2))/2)= &
 595 &                     gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 596                      call metcon(rank,gmet,temp,tmpfac)
 597                      gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
 598 &                     wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 599                    else
 600 !                    ------ Spin-orbit ------
 601                      gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero
 602 !                    Contraction over spins:
 603                      do ispinp=1,nspinortot
 604 !                      => Imaginary part (multiplying by i, then by the Im of amet):
 605                        temp(1,1:((rank+1)*(rank+2))/2)= &
 606 &                       -gxa(2,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
 607                        temp(2,1:((rank+1)*(rank+2))/2)= &
 608 &                       gxa(1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
 609                        amet_lo(:,:)=amet(2,:,:,ispin,ispinp)
 610                        call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
 611                        gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
 612 &                       gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
 613 &                       wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 614 !                      => Real part:
 615                        temp(:,1:((rank+1)*(rank+2))/2)= &
 616 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
 617                        amet_lo(:,:)=amet(1,:,:,ispin,ispinp)
 618                        call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
 619                        gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
 620 &                       gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
 621 &                       wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 622                      end do
 623                    end if
 624 
 625 !                  Eventual tensorial compaction/decompaction for ddk
 626 !                  perturbation:
 627                    if(choice==5 .and. ilang>=2) then
 628                      jjk=1+((ilang-2)*(ilang-1)*ilang)/6
 629                      compact=-1
 630                      temp(:,1:(rank*(rank+1))/2)= &
 631 &                     dgxdt(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj,ispin)
 632                      call ddkten(compact,idir,rank,temp,tmpfac)
 633                      dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)= &
 634 &                     dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)&
 635 &                     +tmpfac(:,1:((rank+1)*(rank+2))/2)
 636                      compact=1
 637                      tmpfac(:,1:((rank+1)*(rank+2))/2)= &
 638 &                     gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)
 639                      call ddkten(compact,idir,rank,temp,tmpfac)
 640                      dgxdtfac(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj)= &
 641 &                     temp(:,1:(rank*(rank+1))/2)
 642                    end if
 643 
 644 !                  Section for strain perturbation
 645 !                  no spin-orbit yet
 646 
 647                    if(choice==3 .and. signs==2) then
 648                      istr=idir
 649                      if(ispinor==1) then
 650 !                      ------ Scalar-Relativistic ------
 651 !                      jjstr is the starting address for dgxds and dgxdsfac
 652                        jjstr=-3+((ilang+1)*(ilang+2)*(ilang+3))/6
 653 !                      diagonal enlk contribution
 654 !                      note sign change (12/05/02)
 655                        if(istr>3) then
 656                          gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero
 657                        else
 658                          gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
 659 &                         -gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)
 660                        end if
 661                        dgxdsfac(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)=zero
 662                        iterm=1
 663                        temp(:,1:((rank+1)*(rank+2))/2)= &
 664 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 665                        call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac)
 666                        dgxdsfac(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)= &
 667 &                       wt(ilang,iproj)*tmpfac(:,1:((rank+3)*(rank+4))/2)
 668                        iterm=2
 669                        temp(:,1:((rank+1)*(rank+2))/2)= &
 670 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 671                        call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac)
 672                        gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
 673 &                       +gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
 674 &                       wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 675                        iterm=3
 676                        temp(:,1:((rank+3)*(rank+4))/2)= &
 677                        dgxds(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)
 678                        call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac)
 679                        gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= &
 680 &                       gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ &
 681 &                       wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 682                      end if
 683 !                    end section for strain perturbation
 684                    end if
 685 
 686 !                  Eventual metric & spinorial contraction from dgxdt to dgxdtfac.
 687 !                  either for the dynamical matrix, or for the application of the
 688 !                  gradient of the operator. The treatment is different for the
 689 !                  scalar-relativistic and spin-orbit parts.
 690                    if ((choice==2.and.signs==2).or. &
 691 &                   (choice==5.and.signs==2).or. &
 692 &                   (choice==4)) then
 693                      mumax=ndgxdtfac;if (choice==5) mumax=1
 694                      do mu=1,mumax
 695                        if(ispinor==1) then
 696 !                        ------ Scalar-Relativistic ------
 697                          temp(:,1:((rank+1)*(rank+2))/2)= &
 698 &                         dgxdt(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 699                          call metcon(rank,gmet,temp,tmpfac)
 700                          dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
 701 &                         wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 702                        else
 703 !                        ------ Spin-orbit ------
 704                          dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero
 705 !                        Contraction over spins:
 706                          do ispinp=1,nspinortot
 707 !                          => Imaginary part (multiplying by i, then by the Im of amet):
 708                            temp(1,1:((rank+1)*(rank+2))/2)= &
 709 &                           -dgxdt(2,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
 710                            temp(2,1:((rank+1)*(rank+2))/2)= &
 711 &                           dgxdt(1,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
 712                            amet_lo(:,:)=amet(2,:,:,ispin,ispinp)
 713                            call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
 714                            dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
 715 &                           dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+&
 716 &                           wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 717 !                          => Real part:
 718                            temp(:,1:((rank+1)*(rank+2))/2)= &
 719 &                           dgxdt(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp)
 720                            amet_lo(:,:)=amet(1,:,:,ispin,ispinp)
 721                            call metcon_so(rank,gmet,amet_lo,temp,tmpfac)
 722                            dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=&
 723 &                           dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+&
 724 &                           wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2)
 725                          end do
 726                        end if
 727                      end do
 728                    end if
 729 
 730 
 731 !                  ----  Accumulate the nonlocal energy.
 732                    do ii=1,ilang2
 733                      jj=ii-1+jjs
 734                      enlk=enlk+(gxafac(1,jj,ia,iproj)*gxa(1,jj,ia,iproj,ispin)&
 735 &                     +gxafac(2,jj,ia,iproj)*gxa(2,jj,ia,iproj,ispin) )
 736                    end do
 737 
 738 !                  ----  Accumulate force contributions if requested.
 739 !                  Note that the contraction of gxa with dgxdt can be done like
 740 !                  a cartesian dot product now because the symmetrically-related
 741 !                  terms are accounted for with weights in gxa.
 742                    if ((choice==2.or.choice==23) .and. signs==1) then
 743                      ishift=0;if (choice==23) ishift=6
 744                      ia5=ia+ia3-1
 745                      do ii=1,ilang2
 746                        jj=ii-1+jjs
 747                        do mu=1,3
 748 !                        (includes also factor of 2 from "2*Re[stuff]")
 749                          indx=mu+3*(ia5-1)+ishift
 750                          enlout(indx)=enlout(indx)+two*&
 751 &                         ( gxafac(1,jj,ia,iproj)*dgxdt(1,mu,jj,ia,iproj,ispin)&
 752 &                         +gxafac(2,jj,ia,iproj)*dgxdt(2,mu,jj,ia,iproj,ispin))
 753                        end do
 754                      end do
 755                    end if
 756 
 757 !                  ----  Accumulate stress tensor contributions if requested.
 758                    if ((choice==3.or.choice==23).and.signs==1) then
 759 !                    1- Compute contractions involving gxa and dgxds:
 760 !                    --- Same formula for Scalar-relativistic and Spin-orbit ---
 761                      if (ilang==1) then
 762                        do ii=1,6
 763                          rank2(ii)=2.d0*&
 764 &                         (gxafac(1,1,ia,iproj)*dgxds(1,ii,ia,iproj,isp)+ &
 765 &                         gxafac(2,1,ia,iproj)*dgxds(2,ii,ia,iproj,isp) )
 766                        end do
 767                      else if (ilang==2) then
 768                        call cont13(gxafac(:,jjs:jjs+2,ia,iproj),&
 769 &                       dgxds(:, 7:16,ia,iproj,isp),rank2)
 770                      else if (ilang==3) then
 771                        call cont24(gxafac(:,jjs:jjs+5,ia,iproj),&
 772 &                       dgxds(:,17:31,ia,iproj,isp),rank2)
 773                      else if (ilang==4) then
 774                        call cont35(gxafac(:,jjs:jjs+9,ia,iproj),&
 775 &                       dgxds(:,32:52,ia,iproj,isp),rank2)
 776                      end if
 777 !                    In all cases add rank2 term into stress tensor
 778                      strsnl(:)=strsnl(:)-rank2(:)
 779 !                    2- Compute contractions involving gxa and gxa:
 780                      if(ispinor==1) then
 781 !                      2a ------ Scalar-Relativistic ------
 782                        if (ilang==2) then
 783                          strsnl(1)=strsnl(1)-wt(ilang,iproj)*2.d0*&
 784 &                         (gxa(1,jjs  ,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispin)+&
 785 &                         gxa(2,jjs  ,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispin))
 786                          strsnl(2)=strsnl(2)-wt(ilang,iproj)*2.d0*&
 787 &                         (gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispin)+&
 788 &                         gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispin))
 789                          strsnl(3)=strsnl(3)-wt(ilang,iproj)*2.d0*&
 790 &                         (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+2,ia,iproj,ispin)+&
 791 &                         gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+2,ia,iproj,ispin))
 792                          strsnl(4)=strsnl(4)-wt(ilang,iproj)*2.d0*&
 793 &                         (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispin)+&
 794 &                         gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispin))
 795                          strsnl(5)=strsnl(5)-wt(ilang,iproj)*2.d0*&
 796 &                         (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispin)+&
 797 &                         gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispin))
 798                          strsnl(6)=strsnl(6)-wt(ilang,iproj)*2.d0*&
 799 &                         (gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispin)+&
 800 &                         gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispin))
 801                        else if (ilang==3) then
 802                          call trace2(gxa(:,jjs:jjs+5,ia,iproj,ispin),gmet,trace)
 803                          call cont22(gxa(:,jjs:jjs+5,ia,iproj,ispin),gmet,rank2)
 804                          do ii=1,6
 805                            strsnl(ii)=strsnl(ii)+wt(ilang,iproj)*&
 806 &                           (2.d0*(trace(1)*gxa(1,jjs+ii-1,ia,iproj,ispin)+&
 807 &                           trace(2)*gxa(2,jjs+ii-1,ia,iproj,ispin))-3.d0*rank2(ii))
 808                          end do
 809                        else if (ilang==4) then
 810                          call cont3(gxa(:,jjs:jjs+9,ia,iproj,ispin),gmet,rank2)
 811                          strsnl(:)=strsnl(:)-wt(ilang,iproj)*rank2(:)
 812                        end if
 813                      else
 814 !                      2b ------ Spin-orbit ------
 815                        do ispinp=1,nspinortot
 816                          if (ilang==3) then
 817                            call cont22so(gxa(:,jjs:jjs+5,ia,iproj,ispin),&
 818 &                           gxa(:,jjs:jjs+5,ia,iproj,ispinp),&
 819 &                           amet(:,:,:,ispin,ispinp),rank2)
 820                            strsnl(:)=strsnl(:)-wt(ilang,iproj)*3.d0*rank2(:)
 821                          else if (ilang==4) then
 822                            call cont33so(gxa(:,jjs:jjs+9,ia,iproj,ispin),&
 823 &                           gxa(:,jjs:jjs+9,ia,iproj,ispinp),&
 824 &                           gmet,amet(:,:,:,ispin,ispinp),rank2)
 825                            strsnl(:)=strsnl(:)-wt(ilang,iproj)*rank2(:)
 826                          end if
 827                        end do
 828                      end if
 829 !                    3- Compute contractions involving gxa and gxa due to
 830 !                    gradients of antisymmetric tensor (amet):
 831 !                    --- Only in case of Spin-orbit ---
 832                      if(ispinor==2) then
 833                        do ispinp=1,nspinortot
 834 !                        Be carefull: no need to compute rank2c(:,1:3) !
 835                          if (ilang==2) then
 836                            rank2c(1,4)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispinp)&
 837 &                           +gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispinp)
 838                            rank2c(2,4)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispinp)&
 839 &                           -gxa(2,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispinp)
 840                            rank2c(1,5)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispinp)&
 841 &                           +gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispinp)
 842                            rank2c(2,5)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispinp)&
 843 &                           -gxa(2,jjs+2,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispinp)
 844                            rank2c(1,6)=gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispinp)&
 845 &                           +gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispinp)
 846                            rank2c(2,6)=gxa(1,jjs+1,ia,iproj,ispin)*gxa(2,jjs  ,ia,iproj,ispinp)&
 847 &                           -gxa(2,jjs+1,ia,iproj,ispin)*gxa(1,jjs  ,ia,iproj,ispinp)
 848                          else if (ilang==3) then
 849                            call cont22cso(gxa(:,jjs:jjs+5,ia,iproj,ispin),&
 850 &                           gxa(:,jjs:jjs+5,ia,iproj,ispinp),&
 851 &                           gmet,rank2c)
 852                          else if (ilang==4) then
 853                            call cont33cso(gxa(:,jjs:jjs+9,ia,iproj,ispin),&
 854 &                           gxa(:,jjs:jjs+9,ia,iproj,ispinp),&
 855 &                           gmet,rank2c)
 856                          end if
 857                          if (ilang>1) then
 858                            do jj=1,3
 859                              do ii=4,6
 860                                strsso(ii,jj)=strsso(ii,jj)-2.d0*wt(ilang,iproj)*&
 861 &                               (pauli(1,ispin,ispinp,jj)*rank2c(2,ii)+&
 862 &                               pauli(2,ispin,ispinp,jj)*rank2c(1,ii))
 863                              end do
 864                            end do
 865                          end if
 866                        end do
 867                      end if
 868 !                    Enf if (choice==3)
 869                    end if
 870 
 871 !                  ----  Accumulate dynamical matrix contributions if requested.
 872                    if (choice==4) then
 873                      ia5=ia+ia3-1
 874                      do ii=1,ilang2
 875                        jj=ii-1+jjs
 876                        do mu=1,6
 877 !                        (includes also factor of 2 from "2*Re[stuff]")
 878                          enlout(mu+6*(ia5-1))=enlout(mu+6*(ia5-1))+two*&
 879 &                         (gxafac(1,jj,ia,iproj)*dgxdt(1,mu+3,jj,ia,iproj,ispin)&
 880 &                         +gxafac(2,jj,ia,iproj)*dgxdt(2,mu+3,jj,ia,iproj,ispin))
 881                        end do
 882                        do mu=1,3
 883                          enlout(mu+6*(ia5-1))=enlout(mu+6*(ia5-1))+two*&
 884 &                         (dgxdtfac(1,mu,jj,ia,iproj)*dgxdt(1,mu,jj,ia,iproj,ispin)&
 885 &                         +dgxdtfac(2,mu,jj,ia,iproj)*dgxdt(2,mu,jj,ia,iproj,ispin))
 886                        end do
 887                        enlout(4+6*(ia5-1))=enlout(4+6*(ia5-1)) +two*&
 888 &                       (dgxdtfac(1,2,jj,ia,iproj)*dgxdt(1,3,jj,ia,iproj,ispin)&
 889 &                       +dgxdtfac(2,2,jj,ia,iproj)*dgxdt(2,3,jj,ia,iproj,ispin))
 890                        enlout(5+6*(ia5-1))=enlout(5+6*(ia5-1)) +two*&
 891 &                       (dgxdtfac(1,1,jj,ia,iproj)*dgxdt(1,3,jj,ia,iproj,ispin)&
 892 &                       +dgxdtfac(2,1,jj,ia,iproj)*dgxdt(2,3,jj,ia,iproj,ispin))
 893                        enlout(6+6*(ia5-1))=enlout(6+6*(ia5-1)) +two*&
 894 &                       (dgxdtfac(1,1,jj,ia,iproj)*dgxdt(1,2,jj,ia,iproj,ispin)&
 895 &                       +dgxdtfac(2,1,jj,ia,iproj)*dgxdt(2,2,jj,ia,iproj,ispin))
 896                      end do
 897                    end if
 898 
 899 !                  ----  Accumulate elastic tensor contributions if requested.
 900 
 901                    if(choice==6) then
 902 !                    XG 081121 : msg to the person who has introduced this CPP option (sorry, I did not have to time to locate who did this ...)
 903 !                    This section of ABINIT should be allowed by default to the user. I have found that on the contrary, the build
 904 !                    system defaults are such that this section is forbidden by default. You might restore this flag if at the same time,
 905 !                    you modify the build system in such a way that by default this section is included, and if the user wants, it can disable it.
 906 !                    #if defined USE_NLSTRAIN_LEGENDRE
 907                      ia5=ia+ia3-1
 908                      jjs1=((ilang)*(ilang+1)*(ilang+2))/6
 909                      jjs2=-3+((ilang+1)*(ilang+2)*(ilang+3))/6
 910                      jjs3=-9+((ilang+2)*(ilang+3)*(ilang+4))/6
 911                      jjs4=-19+((ilang+3)*(ilang+4)*(ilang+5))/6
 912 
 913 !                    Contribution for two diagonal strains (basically, the nonlocal
 914 !                    enlk)
 915                      temp(:,1:((rank+1)*(rank+2))/2)= &
 916 &                     gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 917                      call metcon(rank,gmet,temp,tmpfac)
 918                      e2nldd=zero
 919                      do ii=1,((rank+1)*(rank+2))/2
 920                        e2nldd=e2nldd+&
 921 &                       (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+&
 922 &                       gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii))
 923                      end do
 924 
 925 !                    Terms involving one ucvol derivative (diagonal strain only)
 926 !                    and one derivative of nonlocal operator
 927 !                    Loop over strain index
 928                      do istr2=1,6
 929 
 930 !                      rank->rank+2
 931                        iterm=1
 932                        temp(:,1:((rank+1)*(rank+2))/2)= &
 933 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 934                        call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac)
 935                        e2nl_tmp(istr2)=0.d0
 936                        do ii=1,((rank+3)*(rank+4))/2
 937                          e2nl_tmp(istr2)=e2nl_tmp(istr2)-2.d0*&
 938 &                         (dgxds(1,jjs2-1+ii,ia,iproj,isp)*tmpfac(1,ii)+&
 939 &                         dgxds(2,jjs2-1+ii,ia,iproj,isp)*tmpfac(2,ii))
 940                        end do
 941 !                      rank->rank
 942 
 943                        iterm=2
 944                        temp(:,1:((rank+1)*(rank+2))/2)= &
 945 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)
 946                        call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac)
 947                        do ii=1,((rank+1)*(rank+2))/2
 948                          e2nl_tmp(istr2)=e2nl_tmp(istr2)-&
 949 &                         (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+&
 950 &                         gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii))
 951                        end do
 952 !                      DEBUG
 953 !                      This and subsequent similar debug sections evaluate the
 954 !                      hermitial conjugate of the contraction immeditely above
 955 !                      and the comparison was useful for development purposes.
 956 !                      rank+2->rank
 957 !                      iterm=3
 958 !                      temp(:,1:((rank+3)*(rank+4))/2)= &
 959 !                      dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin)
 960 !                      call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac)
 961 !                      e2nl_tmp(istr2)=0.d0
 962 !                      do ii=1,((rank+1)*(rank+2))/2
 963 !                      e2nl_tmp(istr2)=e2nl_tmp(istr2)-&
 964 !                      &             (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+&
 965 !                      &              gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii))
 966 !                      end do
 967 !                      ENDDEBUG
 968                      end do
 969 
 970 !                    Terms involving two derivatives of the nonlocal operator
 971 !                    Loop over 2nd strain index
 972                      do istr2=1,6
 973 !                      Loop over 1st strain index, upper triangle only
 974                        do istr1=1,6
 975                          iest=istr1+(3*natom+6)*(istr2-1)
 976 
 977 !                        Accumulate terms corresponding to two derivatives of ucvol
 978 !                        (simply the nonlocal energy contributin) for both indices
 979 !                        corresponding to diagonal strains
 980 
 981                          if(istr1<=3 .and. istr2<=3) then
 982                            enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nldd
 983                          end if
 984 
 985 !                        Accumulate terms computed above from 1st derivatives
 986 !                        when one or more indices corresponds to diagonal strain
 987                          if(istr2<=3) then
 988                            enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl_tmp(istr1)
 989                          end if
 990                          if(istr1<=3) then
 991                            enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl_tmp(istr2)
 992                          end if
 993 
 994 !                        rank->rank+4
 995                          call contstr21(istr2,istr1,rank,gmet,gprimd,e2nl,&
 996 &                         gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
 997 &                         d2gxds2(:,jjs4:jjs4-1+((rank+5)*(rank+6))/2,ia,iproj,isp))
 998                          enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
 999 
1000 !                        DEBUG
1001 !                        rank+4->rank
1002 !                        call contstr22(istr2,istr1,rank,gmet,gprimd,e2nl,&
1003 !                        &             d2gxds2(:,jjs4:jjs4-1+((rank+5)*(rank+6))/2,ia,iproj,ispin),&
1004 !                        &             gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1005 !                        enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1006 !                        ENDDEBUG
1007 
1008 !                        rank->rank+2
1009                          call contstr23(istr2,istr1,rank,gmet,gprimd,e2nl,&
1010 &                         gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1011 &                         dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1012                          enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1013 !                        DEBUG
1014 
1015 !                        rank+2->rank
1016 !                        call contstr24(istr2,istr1,rank,gmet,gprimd,e2nl,&
1017 !                        &             dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin),&
1018 !                        &             gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1019 !                        enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1020 !                        ENDDEBUG
1021 
1022 !                        rank+2->rank+2
1023                          if(rank<=2) then
1024                            call contstr25(istr2,istr1,rank,gmet,gprimd,e2nl,&
1025 &                           dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp),&
1026 &                           dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1027                          else
1028                            call contstr25a(istr2,istr1,rank,gmet,gprimd,e2nl,&
1029 &                           dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp),&
1030 &                           dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1031                          end if
1032                          enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1033 
1034 !                        rank->rank
1035                          call contstr26(istr2,istr1,rank,gmet,gprimd,e2nl,&
1036 &                         gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1037 &                         gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1038                          enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl
1039 
1040                        end do !istr1
1041 
1042 !                      Contributions to internal strain (one cartesian strain and one
1043 !                      reduced-coordinate atomic displacement derivative).
1044                        iest=7+3*(ia5-1)+(3*natom+6)*(istr2-1)
1045 
1046 !                      rank->rank+3
1047                        call contistr03(istr2,rank,gmet,gprimd,eisnl,&
1048 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1049 &                       d2gxdis(:,jjs3:jjs3-1+((rank+4)*(rank+5))/2,ia,iproj,isp))
1050                        enlout(iest:iest+2)= enlout(iest:iest+2)&
1051 &                       +wt(ilang,iproj)*eisnl(1:3)
1052 
1053 !                      DEBUG
1054 !                      rank+3->rank
1055 !                      call contistr30(istr2,rank,gmet,gprimd,eisnl,&
1056 !                      &            d2gxdis(:,jjs3:jjs3-1+((rank+4)*(rank+5))/2,ia,iproj,ispin),&
1057 !                      &            gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1058 !                      enlout(iest:iest+2)= enlout(iest:iest+2)&
1059 !                      &            +wt(ilang,iproj)*eisnl(1:3)
1060 !                      ENDDEBUG
1061 
1062 !                      rank+1->rank+2
1063                        call contistr12(istr2,rank,gmet,gprimd,eisnl,&
1064 &                       dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,isp),&
1065 &                       dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp))
1066                        enlout(iest:iest+2)= enlout(iest:iest+2)&
1067 &                       +wt(ilang,iproj)*eisnl(1:3)
1068 
1069 !                      DEBUG
1070 !                      rank+2->rank+1
1071 !                      call contistr21(istr2,rank,gmet,gprimd,eisnl,&
1072 !                      &            dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin),&
1073 !                      &            dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,ispin))
1074 !                      enlout(iest:iest+2)= enlout(iest:iest+2)&
1075 !                      &            +wt(ilang,iproj)*eisnl(1:3)
1076 !                      ENDDEBUG
1077 
1078 !                      rank->rank+1
1079                        call contistr01(istr2,rank,gmet,gprimd,eisnl,&
1080 &                       gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),&
1081 &                       dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,isp))
1082                        enlout(iest:iest+2)= enlout(iest:iest+2)&
1083 &                       +wt(ilang,iproj)*eisnl(1:3)
1084 !
1085 !                      DEBUG
1086 !                      rank+1->rank
1087 !                      call contistr10(istr2,rank,gmet,gprimd,eisnl,&
1088 !                      &            dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,ispin),&
1089 !                      &            gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin))
1090 !                      enlout(iest:iest+2)= enlout(iest:iest+2)&
1091 !                      &            +wt(ilang,iproj)*eisnl(1:3)
1092 !                      ENDDEBUG
1093 
1094                      end do !istr2
1095                    end if !choice==6
1096 
1097 !                  End loop over iproj:
1098                  end do
1099 !                End condition of existence of a reference state:
1100                end if
1101 
1102 !              End loop over ilang:
1103              end do
1104 
1105 !            End loop over ia:
1106            end do
1107 
1108 !          Operate with the non-local potential on the projected scalars,
1109 !          in order to get matrix element [NOT needed if only force or stress
1110 !          or dynamical matrix is being computed].
1111 
1112            if (signs==2) then
1113              if(nloalg(2)<=0 .and. choice==2)then
1114 !              Prepare the phase k+q factors for the atoms between ia3 and ia4:
1115 !              they were not prepared previously for nloalg(2)<=0 and choice==2.
1116                call ph1d3d(ia3,ia4,kgout,matblk,natom,npwout,&
1117 &               n1,n2,n3,phkxredout,ph1d,ph3dout)
1118              end if
1119 
1120 !            call timab(74+choice,1,tsec)
1121              sign=-1
1122 !            The duplication of the opernl calls has the purpose to avoid
1123 !            copying vectout/vectout_s
1124              if(ispin==1.or.nspinor==1)then
1125                if( nloalg(1)==2) then
1126                  call opernl2(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1127 &                 ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1128 &                 jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1129 &                 mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1130 &                 ntypat,ph3dout,sign,vectout)
1131                else if( nloalg(1)==3) then
1132                  call opernl3(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1133 &                 ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1134 &                 jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1135 &                 mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1136 &                 ntypat,ph3dout,sign,vectout)
1137                else if( nloalg(1)==4) then
1138                  call opernl4b(choice,dgxdsfac,dgxdtfac,ffnlout,gmet,gxafac,&
1139 &                 ia3,idir,indlmn_s,ispinor,itypat,jproj,kgout,kpgout,kptout,&
1140 &                 lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,&
1141 &                 dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,ntypat,ph3dout,vectout)
1142                end if
1143              else  ! if ispin == 2
1144                vectout_s(:,:)=zero
1145                if( nloalg(1)==2) then
1146                  call opernl2(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1147 &                 ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1148 &                 jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1149 &                 mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1150 &                 ntypat,ph3dout,sign,vectout_s)
1151                else if( nloalg(1)==3) then
1152                  call opernl3(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,&
1153 &                 ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,&
1154 &                 jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1155 &                 mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,&
1156 &                 ntypat,ph3dout,sign,vectout_s)
1157                else if( nloalg(1)==4) then
1158                  call opernl4b(choice,dgxds,dgxdtfac,ffnlout,gmet,gxafac,&
1159 &                 ia3,idir,indlmn_s,ispinor,itypat,jproj,kgout,kpgout,kptout,&
1160 &                 lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,&
1161 &                 dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,ntypat,ph3dout,vectout_s)
1162                end if
1163                vectout(1,1+npwout:2*npwout)=&
1164 &               vectout(1,1+npwout:2*npwout)+vectout_s(1,1:npwout)
1165                vectout(2,1+npwout:2*npwout)=&
1166 &               vectout(2,1+npwout:2*npwout)+vectout_s(2,1:npwout)
1167 
1168              end if ! end ispin if
1169            end if ! end signs==2 if
1170 
1171 !          End loops over spins:
1172          end do
1173 
1174 !        End condition of existence of a non-local part for that type of atom:
1175        end if
1176 
1177 !      End loop over ispinor:
1178      end do
1179 
1180 !    End sum on atom subset loop, over ia3:
1181    end do
1182 
1183 !  End atom type loop, over itypat:
1184    ia1=ia2+1
1185  end do
1186 
1187 !De-allocate temporary space.
1188  ABI_FREE(ekb_s)
1189  ABI_FREE(gxa)
1190  ABI_FREE(gxafac)
1191  ABI_FREE(dgxds)
1192  ABI_FREE(dgxdt)
1193  ABI_FREE(dgxdtfac)
1194  ABI_FREE(wt)
1195  ABI_FREE(jproj)
1196  ABI_FREE(temp)
1197  ABI_FREE(tmpfac)
1198  ABI_FREE(dgxdsfac)
1199  ABI_FREE(indlmn_s)
1200  !if(choice==6)  then
1201  ABI_FREE(dgxdis)
1202  ABI_FREE(d2gxdis)
1203  ABI_FREE(d2gxds2)
1204  !end if
1205  !if(nspinor==2) then
1206  ABI_FREE(dgxds_s)
1207  ABI_FREE(dgxdt_s)
1208  ABI_FREE(gxa_s)
1209  !end if
1210  !if(nspinor==2.and.choice==6) then
1211  ABI_FREE(dgxdis_s)
1212  ABI_FREE(d2gxdis_s)
1213  ABI_FREE(d2gxds2_s)
1214  !end if
1215  if (mpssoang>mpsang)  then
1216    ABI_FREE(pauli)
1217  end if
1218 
1219 !Restore the original content of the vectin array.
1220 !Note that only the first part was modified
1221  if(istwf_k/=1) then
1222    call scalewf_nonlop(istwf_k,mpi_enreg,npwin,2,vectin)
1223  end if
1224 
1225  if (nspinor==2)  then
1226    ABI_FREE(vectin_s)
1227    ABI_FREE(vectout_s)
1228  end if
1229 
1230  if (mpi_enreg%paral_spinor==1) then
1231    if (size(enlout)>0) call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr)
1232    call xmpi_sum(strsnl,mpi_enreg%comm_spinor,ierr)
1233    call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr)
1234    call xmpi_sum(strsso,mpi_enreg%comm_spinor,ierr)
1235  end if
1236 
1237 !Save non-local energy
1238  if((choice==1).and.size(enlout)>0) enlout(1)=enlk ! on test v4/93 size(enlout) can be zero (PMA)
1239 
1240 !Do final manipulations to produce strain gradients for
1241 !stress tensor, in cartesian coordinates
1242  if ((choice==3.or.choice==23) .and. signs==1) then
1243 !  Convert strsnl from reduced to cartesian coordinates
1244    strsnl_out(:)=0.d0
1245    call strconv(strsnl,gprimd,strsnl_out)
1246    strsnl(:) = strsnl_out(:)
1247 !  Add diagonal part (fill up first 6 components of enlout with
1248 !  these gradients; elements 7,8,9 of enlout are not used)
1249    enlout(1)=strsnl(1)-enlk
1250    enlout(2)=strsnl(2)-enlk
1251    enlout(3)=strsnl(3)-enlk
1252    enlout(4)=strsnl(4)
1253    enlout(5)=strsnl(5)
1254    enlout(6)=strsnl(6)
1255 !  Eventually, add spin-orbit part due to gradients of
1256 !  antisymmetric tensor
1257    if (mpssoang>mpsang) then
1258      call strsocv(strsso,gprimd,strssoc)
1259      enlout(1:6)=enlout(1:6)+strssoc(:)
1260    end if
1261  end if
1262 
1263 !DEBUG
1264 !write(std_out,*)' nonlop_pl: exit '
1265 !ENDDEBUG
1266 
1267 contains

ABINIT/scalewf_nonlop [ Functions ]

[ Top ] [ Functions ]

NAME

 scalewf_nonlop

FUNCTION

 At the start of nonlop (or similar routines), as well as its end,
 the wavefunctions, when stored with istwfk/=2,
 need to be scaled (by a factor of 2 or 1/2),
 except for the G=0 component.
 Only the first spinor component is to be modified.

INPUTS

  istwf_k=storage mode of the vector
  mpi_enreg=information about MPI parallelization
  npw=number of planewaves
  option=1 multiply by 2
        =2 multiply by 1/2

OUTPUT

  (see side effects)

SIDE EFFECTS

  vect(2,npw)=vector that is rescaled

NOTES

  XG030513 : MPIWF One should pay attention to the
  G=0 component, that will be only one one proc...

SOURCE

1431 subroutine scalewf_nonlop(istwf_k,mpi_enreg,npw,option,vect)
1432 
1433 !Arguments ------------------------------------
1434 !This type is defined in defs_mpi
1435 !scalars
1436  integer,intent(in) :: istwf_k,npw,option
1437  type(MPI_type),intent(in) :: mpi_enreg
1438 !arrays
1439  real(dp),intent(inout) :: vect(2,npw)
1440 
1441 !Local variables-------------------------------
1442 !scalars
1443  integer :: ipw
1444  real(dp) :: scale
1445  character(len=500) :: msg
1446 
1447 ! *************************************************************************
1448 
1449  DBG_ENTER("COLL")
1450 
1451  if(istwf_k/=1)then
1452 
1453    if(option/=1 .and. option/=2)then
1454      write(msg,'(a,a,a,i0)')&
1455      'The argument option should be 1 or 2,',ch10,&
1456      'however, option=',option
1457      ABI_BUG(msg)
1458    end if
1459 
1460    scale=two
1461    if(option==2)scale=half
1462 
1463 !  Storage for the Gamma point. The component of the G=0 vector
1464 !  should not be scaled, and no G=0 imaginary part is allowed.
1465    if(istwf_k==2)then
1466      if (mpi_enreg%me_g0==1) then
1467        vect(2,1)=zero
1468 !$OMP PARALLEL DO
1469        do ipw=2,npw
1470          vect(1,ipw)=scale*vect(1,ipw)
1471          vect(2,ipw)=scale*vect(2,ipw)
1472        end do
1473 !$OMP END PARALLEL DO
1474      else
1475 !$OMP PARALLEL DO
1476        do ipw=1,npw
1477          vect(1,ipw)=scale*vect(1,ipw)
1478          vect(2,ipw)=scale*vect(2,ipw)
1479        end do
1480 !$OMP END PARALLEL DO
1481      end if
1482    end if
1483 
1484 !  Other storage modes, for k points with time-reversal symmetry.
1485 !  All components should be scaled.
1486    if(istwf_k>2)then
1487 !$OMP PARALLEL DO
1488      do ipw=1,npw
1489        vect(1,ipw)=scale*vect(1,ipw)
1490        vect(2,ipw)=scale*vect(2,ipw)
1491      end do
1492 !$OMP END PARALLEL DO
1493    end if
1494 
1495  end if ! istwf_k/=1
1496 
1497  DBG_EXIT("COLL")
1498 
1499 end subroutine scalewf_nonlop

ABINIT/strsocv [ Functions ]

[ Top ] [ Functions ]

NAME

 strsocv

FUNCTION

 Convert from antisymmetric storage mode 3x3x3 rank3 tensor in reduced
 coordinates "red" to symmetric storage mode 3x3 rank2 tensor in
 cartesian coordinates "cart", using metric tensor "gprimd".

INPUTS

  red(6,3)=3x3x3 tensor in antisymmetric storage mode,
           reduced coordinates
  gprimd(3,3)=reciprocal space dimensional primitive translations

OUTPUT

  cart(6)=3x3 tensor in symmetric storage mode,
          cartesian coordinates

NOTES

 This routine is used to compute spin-orbit stress tensor.

 red is antisymmetric in first two indices and stored
     as 11 22 33 32 31 21.
 cart is symmetric and stored as 11 22 33 32 31 21.

 cart(1,1) & = &        & red(i,j,2) G(3,i) G(1,j) + red(i,j,3) G(1,i) G(2,j) \nonumber
 cart(2,2) & = &        & red(i,j,1) G(2,i) G(3,j) + red(i,j,3) G(1,i) G(2,j) \nonumber
 cart(3,3) & = &        & red(i,j,1) G(2,i) G(3,j) + red(i,j,2) G(3,i) G(1,j) \nonumber
 cart(3,2) & = &  0.5 ( & red(i,j,3) G(1,i) G(3,j) + red(i,j,2) G(2,i) G(1,j)) \nonumber
 cart(3,1) & = &  0.5 ( & red(i,j,3) G(3,i) G(2,j) + red(i,j,1) G(2,i) G(1,j)) \nonumber
 cart(2,1) & = &  0.5 ( & red(i,j,2) G(3,i) G(2,j) + red(i,j,1) G(1,i) G(3,j))
 \end{eqnarray} }}

SOURCE

1353 subroutine strsocv(red,gprimd,cart)
1354 
1355 !Arguments ------------------------------------
1356 !arrays
1357  real(dp),intent(in) :: gprimd(3,3),red(6,3)
1358  real(dp),intent(out) :: cart(6)
1359 
1360 !Local variables-------------------------------
1361 !scalars
1362  integer :: ii,jj
1363 !arrays
1364  real(dp) :: work(3,3,3)
1365 
1366 ! *************************************************************************
1367 
1368  do ii=1,3
1369    work(1,1,ii)=0.d0
1370    work(2,2,ii)=0.d0
1371    work(3,3,ii)=0.d0
1372    work(3,2,ii)=red(4,ii) ; work(2,3,ii)=-red(4,ii)
1373    work(3,1,ii)=red(5,ii) ; work(1,3,ii)=-red(5,ii)
1374    work(2,1,ii)=red(6,ii) ; work(1,2,ii)=-red(6,ii)
1375  end do
1376 
1377  cart(:)=0.d0
1378  do jj=1,3
1379    do ii=1,3
1380      cart(1)=cart(1)+work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)&
1381 &     +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj)
1382      cart(2)=cart(2)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)&
1383 &     +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj)
1384      cart(3)=cart(3)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)&
1385 &     +work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)
1386      cart(4)=cart(4)+work(ii,jj,3)*gprimd(1,ii)*gprimd(3,jj)&
1387 &     +work(ii,jj,2)*gprimd(2,ii)*gprimd(1,jj)
1388      cart(5)=cart(5)+work(ii,jj,3)*gprimd(3,ii)*gprimd(2,jj)&
1389 &     +work(ii,jj,1)*gprimd(2,ii)*gprimd(1,jj)
1390      cart(6)=cart(6)+work(ii,jj,2)*gprimd(3,ii)*gprimd(2,jj)&
1391 &     +work(ii,jj,1)*gprimd(1,ii)*gprimd(3,jj)
1392    end do
1393  end do
1394  cart(4)=0.5d0*cart(4)
1395  cart(5)=0.5d0*cart(5)
1396  cart(6)=0.5d0*cart(6)
1397 
1398 end subroutine strsocv

ABINIT/trace2 [ Functions ]

[ Top ] [ Functions ]

NAME

 trace2

FUNCTION

 Sum indices to compute trace of rank 2 tensor gxa related to l=2
 $trace=sum_{i,j} {gxa(i,j) gmet(i,j)}$

INPUTS

  gxa(2,6) = $sum_{G} {e^(2 \pi i G cdot t} {{f_2(|k+G|)} \over {|k+G|^2}} (k+G) cdot (k+G) C(G_{nk})}$
  gmet(3,3)=(symmetric) metric tensor in reciprocal space (bohr^-2)

OUTPUT

  trace(2)=sum_{i,j} {gxa(i,j) gmet(i,j)}$ (Re and Im).

NOTES

 Here index 6 refers to vector components
 of (k+G) but note tensor is symmetric=>only 6 components.
 The components are given in the order 11 22 33 32 31 21.
 The initial 2 handles the Re and Im parts.

SOURCE

1293 subroutine trace2(gxa,gmet,trace)
1294 
1295 !Arguments ------------------------------------
1296 !arrays
1297  real(dp),intent(in) :: gmet(3,3),gxa(2,6)
1298  real(dp),intent(out) :: trace(2)
1299 
1300 !Local variables-------------------------------
1301 !scalars
1302  integer :: reim
1303 
1304 ! *************************************************************************
1305 
1306 !Write out index summation, Re and Im parts
1307  do reim=1,2
1308    trace(reim)=gxa(reim,1)*gmet(1,1)+gxa(reim,2)*gmet(2,2)+&
1309 &   gxa(reim,3)*gmet(3,3)+&
1310 &   2.0d0*(gxa(reim,4)*gmet(3,2)+gxa(reim,5)*gmet(3,1)+&
1311 &   gxa(reim,6)*gmet(2,1))
1312  end do
1313 
1314 end subroutine trace2