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