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.

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1607 subroutine ddkten(compact,idir,rank,temp,tmpfac)
1608 
1609 
1610 !This section has been created automatically by the script Abilint (TD).
1611 !Do not modify the following lines by hand.
1612 #undef ABI_FUNC
1613 #define ABI_FUNC 'ddkten'
1614 !End of the abilint section
1615 
1616  implicit none
1617 
1618 !Arguments ------------------------------------
1619 !scalars
1620  integer,intent(in) :: compact,idir,rank
1621 !arrays
1622  real(dp),intent(inout) :: temp(2,(rank*(rank+1))/2)
1623  real(dp),intent(inout) :: tmpfac(2,((rank+1)*(rank+2))/2)
1624 
1625 !Local variables-------------------------------
1626 !scalars
1627  character(len=500) :: message
1628 
1629 ! *************************************************************************
1630 
1631  if(rank/=1 .and. rank/=2 .and. rank/=3)then
1632    write(message, '(a,i10,a,a,a)' )&
1633 &   'Input rank=',rank,' not allowed.',ch10,&
1634 &   'Possible values are 1,2,3 only.'
1635    MSG_BUG(message)
1636  end if
1637 
1638 !Take care of p angular momentum
1639  if(rank==1)then
1640 
1641 !  Compaction tmpfac -> temp
1642    if(compact==1)then
1643      temp(:,1)=tmpfac(:,idir)
1644 
1645 !    Decompaction temp -> tmpfac
1646    else
1647      tmpfac(:,1:3)=0.0d0
1648      tmpfac(:,idir)=temp(:,1)
1649    end if
1650 
1651 !  Take care of d angular momentum
1652 !  rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1653 
1654  else if(rank==2)then
1655 
1656 !  Compaction tmpfac -> temp
1657    if(compact==1)then
1658      if(idir==1)then
1659 !      Count the number of non-zero derivatives with respect to k(idir)
1660 !      The factor of 2 on the diagonal comes from the derivative with
1661 !      respect to the first K then to the second K
1662        temp(:,1)=2.0d0*tmpfac(:,1); temp(:,2)=tmpfac(:,6); temp(:,3)=tmpfac(:,5)
1663      else if(idir==2)then
1664        temp(:,2)=2.0d0*tmpfac(:,2); temp(:,1)=tmpfac(:,6); temp(:,3)=tmpfac(:,4)
1665      else if(idir==3)then
1666        temp(:,3)=2.0d0*tmpfac(:,3); temp(:,1)=tmpfac(:,5); temp(:,2)=tmpfac(:,4)
1667      end if
1668 !    Decompaction temp -> tmpfac
1669    else
1670      tmpfac(:,1:6)=0.0d0
1671      tmpfac(:,idir)=2.0d0*temp(:,idir)
1672      if(idir==1)then
1673        tmpfac(:,5)=temp(:,3); tmpfac(:,6)=temp(:,2)
1674      else if(idir==2)then
1675        tmpfac(:,4)=temp(:,3); tmpfac(:,6)=temp(:,1)
1676      else if(idir==3)then
1677        tmpfac(:,4)=temp(:,2); tmpfac(:,5)=temp(:,1)
1678      end if
1679    end if
1680 
1681 !  Take care of f angular momentum
1682  else if(rank==3)then
1683 !  rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10
1684 !  rank=2 11->1 22->2 33->3 32->4 31->5 21->6
1685 
1686 !  Compaction tmpfac -> temp
1687    if(compact==1)then
1688      if(idir==1)then
1689 !      Count the number of non-zero derivatives with respect to k(idir)
1690        temp(:,1)=3.0d0*tmpfac(:,1)
1691        temp(:,2:4)=tmpfac(:,2:4)
1692        temp(:,5:6)=2.0d0*tmpfac(:,5:6)
1693      else if(idir==2)then
1694        temp(:,6)=2.0d0*tmpfac(:,2)
1695        temp(:,4)=2.0d0*tmpfac(:,9)
1696        temp(:,5)=tmpfac(:,4)
1697        temp(:,1)=tmpfac(:,6)
1698        temp(:,3)=tmpfac(:,8)
1699        temp(:,2)=3.0d0*tmpfac(:,7)
1700      else if(idir==3)then
1701        temp(:,3)=3.0d0*tmpfac(:,10)
1702        temp(:,5)=2.0d0*tmpfac(:,3)
1703        temp(:,4)=2.0d0*tmpfac(:,8)
1704        temp(:,6)=tmpfac(:,4)
1705        temp(:,1)=tmpfac(:,5)
1706        temp(:,2)=tmpfac(:,9)
1707      end if
1708 !    Decompaction temp -> tmpfac
1709    else
1710      tmpfac(:,1:10)=0.0d0
1711      if(idir==1)then
1712        tmpfac(:,1)=3.0d0*temp(:,1)
1713        tmpfac(:,2:4)=temp(:,2:4)
1714        tmpfac(:,5:6)=2.0d0*temp(:,5:6)
1715      else if(idir==2)then
1716        tmpfac(:,2)=2.0d0*temp(:,6)
1717        tmpfac(:,9)=2.0d0*temp(:,4)
1718        tmpfac(:,4)=temp(:,5)
1719        tmpfac(:,6)=temp(:,1)
1720        tmpfac(:,8)=temp(:,3)
1721        tmpfac(:,7)=3.0d0*temp(:,2)
1722      else if(idir==3)then
1723        tmpfac(:,10)=3.0d0*temp(:,3)
1724        tmpfac(:,3)=2.0d0*temp(:,5)
1725        tmpfac(:,8)=2.0d0*temp(:,4)
1726        tmpfac(:,4)=temp(:,6)
1727        tmpfac(:,5)=temp(:,1)
1728        tmpfac(:,9)=temp(:,2)
1729      end if
1730    end if
1731 
1732  end if
1733 
1734 end subroutine ddkten

ABINIT/m_nonlop_pl [ Modules ]

[ Top ] [ Modules ]

NAME

  nonlop_pl

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_nonlop_pl
27 
28  use defs_basis
29  use defs_abitypes
30  use m_errors
31  use m_abicore
32  use m_xmpi
33  use m_contistr01
34  use m_contistr03
35  use m_contistr12
36  use m_contstr21
37  use m_contstr23
38  use m_contstr25
39  use m_contstr25a
40  use m_contstr26
41  use m_metstr
42  use m_opernl
43 
44  use m_geometry,   only : strconv
45  use m_kg,         only : ph1d3d
46  use m_contract,   only : cont22cso, cont22so, cont24, cont33cso, cont33so, cont35, cont22, cont3, cont13, &
47                           metcon, metcon_so, metric_so
48  implicit none
49 
50  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

PARENTS

      nonlop

CHILDREN

      cont13,cont22,cont22cso,cont22so,cont24,cont3,cont33cso,cont33so,cont35
      contistr01,contistr03,contistr12,contstr21,contstr23,contstr25
      contstr25a,contstr26,ddkten,metcon,metcon_so,metric_so,metstr,opernl2
      opernl3,opernl4a,opernl4b,ph1d3d,scalewf_nonlop,strconv,strsocv,trace2
      xmpi_sum

SOURCE

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

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1488 subroutine scalewf_nonlop(istwf_k,mpi_enreg,npw,option,vect)
1489 
1490 
1491 !This section has been created automatically by the script Abilint (TD).
1492 !Do not modify the following lines by hand.
1493 #undef ABI_FUNC
1494 #define ABI_FUNC 'scalewf_nonlop'
1495 !End of the abilint section
1496 
1497  implicit none
1498 
1499 !Arguments ------------------------------------
1500 !This type is defined in defs_mpi
1501 !scalars
1502  integer,intent(in) :: istwf_k,npw,option
1503  type(MPI_type),intent(in) :: mpi_enreg
1504 !arrays
1505  real(dp),intent(inout) :: vect(2,npw)
1506 
1507 !Local variables-------------------------------
1508 !scalars
1509  integer :: ipw
1510  real(dp) :: scale
1511  character(len=500) :: message
1512 
1513 ! *************************************************************************
1514 
1515  DBG_ENTER("COLL")
1516 
1517  if(istwf_k/=1)then
1518 
1519    if(option/=1 .and. option/=2)then
1520      write(message,'(a,a,a,i0)')&
1521 &     'The argument option should be 1 or 2,',ch10,&
1522 &     'however, option=',option
1523      MSG_BUG(message)
1524    end if
1525 
1526    scale=two
1527    if(option==2)scale=half
1528 
1529 !  Storage for the Gamma point. The component of the G=0 vector
1530 !  should not be scaled, and no G=0 imaginary part is allowed.
1531    if(istwf_k==2)then
1532      if (mpi_enreg%me_g0==1) then
1533        vect(2,1)=zero
1534 !$OMP PARALLEL DO
1535        do ipw=2,npw
1536          vect(1,ipw)=scale*vect(1,ipw)
1537          vect(2,ipw)=scale*vect(2,ipw)
1538        end do
1539 !$OMP END PARALLEL DO
1540      else
1541 !$OMP PARALLEL DO
1542        do ipw=1,npw
1543          vect(1,ipw)=scale*vect(1,ipw)
1544          vect(2,ipw)=scale*vect(2,ipw)
1545        end do
1546 !$OMP END PARALLEL DO
1547      end if
1548    end if
1549 
1550 !  Other storage modes, for k points with time-reversal symmetry.
1551 !  All components should be scaled.
1552    if(istwf_k>2)then
1553 !$OMP PARALLEL DO
1554      do ipw=1,npw
1555        vect(1,ipw)=scale*vect(1,ipw)
1556        vect(2,ipw)=scale*vect(2,ipw)
1557      end do
1558 !$OMP END PARALLEL DO
1559    end if
1560 
1561  end if ! istwf_k/=1
1562 
1563  DBG_EXIT("COLL")
1564 
1565 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} }}

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1396 subroutine strsocv(red,gprimd,cart)
1397 
1398 
1399 !This section has been created automatically by the script Abilint (TD).
1400 !Do not modify the following lines by hand.
1401 #undef ABI_FUNC
1402 #define ABI_FUNC 'strsocv'
1403 !End of the abilint section
1404 
1405  implicit none
1406 
1407 !Arguments ------------------------------------
1408 !arrays
1409  real(dp),intent(in) :: gprimd(3,3),red(6,3)
1410  real(dp),intent(out) :: cart(6)
1411 
1412 !Local variables-------------------------------
1413 !scalars
1414  integer :: ii,jj
1415 !arrays
1416  real(dp) :: work(3,3,3)
1417 
1418 ! *************************************************************************
1419 
1420  do ii=1,3
1421    work(1,1,ii)=0.d0
1422    work(2,2,ii)=0.d0
1423    work(3,3,ii)=0.d0
1424    work(3,2,ii)=red(4,ii) ; work(2,3,ii)=-red(4,ii)
1425    work(3,1,ii)=red(5,ii) ; work(1,3,ii)=-red(5,ii)
1426    work(2,1,ii)=red(6,ii) ; work(1,2,ii)=-red(6,ii)
1427  end do
1428 
1429  cart(:)=0.d0
1430  do jj=1,3
1431    do ii=1,3
1432      cart(1)=cart(1)+work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)&
1433 &     +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj)
1434      cart(2)=cart(2)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)&
1435 &     +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj)
1436      cart(3)=cart(3)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)&
1437 &     +work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)
1438      cart(4)=cart(4)+work(ii,jj,3)*gprimd(1,ii)*gprimd(3,jj)&
1439 &     +work(ii,jj,2)*gprimd(2,ii)*gprimd(1,jj)
1440      cart(5)=cart(5)+work(ii,jj,3)*gprimd(3,ii)*gprimd(2,jj)&
1441 &     +work(ii,jj,1)*gprimd(2,ii)*gprimd(1,jj)
1442      cart(6)=cart(6)+work(ii,jj,2)*gprimd(3,ii)*gprimd(2,jj)&
1443 &     +work(ii,jj,1)*gprimd(1,ii)*gprimd(3,jj)
1444    end do
1445  end do
1446  cart(4)=0.5d0*cart(4)
1447  cart(5)=0.5d0*cart(5)
1448  cart(6)=0.5d0*cart(6)
1449 
1450 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.

PARENTS

      nonlop_pl

CHILDREN

SOURCE

1322 subroutine trace2(gxa,gmet,trace)
1323 
1324 
1325 !This section has been created automatically by the script Abilint (TD).
1326 !Do not modify the following lines by hand.
1327 #undef ABI_FUNC
1328 #define ABI_FUNC 'trace2'
1329 !End of the abilint section
1330 
1331  implicit none
1332 
1333 !Arguments ------------------------------------
1334 !arrays
1335  real(dp),intent(in) :: gmet(3,3),gxa(2,6)
1336  real(dp),intent(out) :: trace(2)
1337 
1338 !Local variables-------------------------------
1339 !scalars
1340  integer :: reim
1341 
1342 ! *************************************************************************
1343 
1344 !Write out index summation, Re and Im parts
1345  do reim=1,2
1346    trace(reim)=gxa(reim,1)*gmet(1,1)+gxa(reim,2)*gmet(2,2)+&
1347 &   gxa(reim,3)*gmet(3,3)+&
1348 &   2.0d0*(gxa(reim,4)*gmet(3,2)+gxa(reim,5)*gmet(3,1)+&
1349 &   gxa(reim,6)*gmet(2,1))
1350  end do
1351 
1352 end subroutine trace2