TABLE OF CONTENTS
- ABINIT/ddkten
- ABINIT/m_nonlop_pl
- ABINIT/nonlop_pl
- ABINIT/scalewf_nonlop
- ABINIT/strsocv
- ABINIT/trace2
ABINIT/ddkten [ Functions ]
NAME
ddkten
FUNCTION
Compact or decompact the tensors related to the ffnl(:,1,...) part of the ddk operator, taking into account the direction of the ddk perturbation.
INPUTS
compact= if 1, compact from tmpfac idir=direction of the ddk perturbation rank=0,1,2, or 3 = rank of tmpfac tensor, also angular momentum (=l)
OUTPUT
(see side effects)
SIDE EFFECTS
Input/Output: temp(2,(rank*(rank+1))/2)=compacted tensor for l=1, just a scalar for l=2, a vector tmpfac(2,(rank+1)*(rank+2)/2)=decompacted tensor for l=1, a vector for l=2, a symmetric matrix, stored as (1 . .) (6 2 .) (5 4 3)
NOTES
For l=0, there is no contribution.
SOURCE
1536 subroutine ddkten(compact,idir,rank,temp,tmpfac) 1537 1538 !Arguments ------------------------------------ 1539 !scalars 1540 integer,intent(in) :: compact,idir,rank 1541 !arrays 1542 real(dp),intent(inout) :: temp(2,(rank*(rank+1))/2) 1543 real(dp),intent(inout) :: tmpfac(2,((rank+1)*(rank+2))/2) 1544 1545 !Local variables------------------------------- 1546 !scalars 1547 character(len=500) :: msg 1548 1549 ! ************************************************************************* 1550 1551 if(rank/=1 .and. rank/=2 .and. rank/=3)then 1552 write(msg, '(a,i10,a,a,a)' )& 1553 'Input rank=',rank,' not allowed.',ch10,& 1554 'Possible values are 1,2,3 only.' 1555 ABI_BUG(msg) 1556 end if 1557 1558 !Take care of p angular momentum 1559 if(rank==1)then 1560 1561 ! Compaction tmpfac -> temp 1562 if(compact==1)then 1563 temp(:,1)=tmpfac(:,idir) 1564 1565 ! Decompaction temp -> tmpfac 1566 else 1567 tmpfac(:,1:3)=0.0d0 1568 tmpfac(:,idir)=temp(:,1) 1569 end if 1570 1571 ! Take care of d angular momentum 1572 ! rank=2 11->1 22->2 33->3 32->4 31->5 21->6 1573 1574 else if(rank==2)then 1575 1576 ! Compaction tmpfac -> temp 1577 if(compact==1)then 1578 if(idir==1)then 1579 ! Count the number of non-zero derivatives with respect to k(idir) 1580 ! The factor of 2 on the diagonal comes from the derivative with 1581 ! respect to the first K then to the second K 1582 temp(:,1)=2.0d0*tmpfac(:,1); temp(:,2)=tmpfac(:,6); temp(:,3)=tmpfac(:,5) 1583 else if(idir==2)then 1584 temp(:,2)=2.0d0*tmpfac(:,2); temp(:,1)=tmpfac(:,6); temp(:,3)=tmpfac(:,4) 1585 else if(idir==3)then 1586 temp(:,3)=2.0d0*tmpfac(:,3); temp(:,1)=tmpfac(:,5); temp(:,2)=tmpfac(:,4) 1587 end if 1588 ! Decompaction temp -> tmpfac 1589 else 1590 tmpfac(:,1:6)=0.0d0 1591 tmpfac(:,idir)=2.0d0*temp(:,idir) 1592 if(idir==1)then 1593 tmpfac(:,5)=temp(:,3); tmpfac(:,6)=temp(:,2) 1594 else if(idir==2)then 1595 tmpfac(:,4)=temp(:,3); tmpfac(:,6)=temp(:,1) 1596 else if(idir==3)then 1597 tmpfac(:,4)=temp(:,2); tmpfac(:,5)=temp(:,1) 1598 end if 1599 end if 1600 1601 ! Take care of f angular momentum 1602 else if(rank==3)then 1603 ! rank=3 111->1 221->2 331->3 321->4 311->5 211->6 222->7 332->8 322->9 333->10 1604 ! rank=2 11->1 22->2 33->3 32->4 31->5 21->6 1605 1606 ! Compaction tmpfac -> temp 1607 if(compact==1)then 1608 if(idir==1)then 1609 ! Count the number of non-zero derivatives with respect to k(idir) 1610 temp(:,1)=3.0d0*tmpfac(:,1) 1611 temp(:,2:4)=tmpfac(:,2:4) 1612 temp(:,5:6)=2.0d0*tmpfac(:,5:6) 1613 else if(idir==2)then 1614 temp(:,6)=2.0d0*tmpfac(:,2) 1615 temp(:,4)=2.0d0*tmpfac(:,9) 1616 temp(:,5)=tmpfac(:,4) 1617 temp(:,1)=tmpfac(:,6) 1618 temp(:,3)=tmpfac(:,8) 1619 temp(:,2)=3.0d0*tmpfac(:,7) 1620 else if(idir==3)then 1621 temp(:,3)=3.0d0*tmpfac(:,10) 1622 temp(:,5)=2.0d0*tmpfac(:,3) 1623 temp(:,4)=2.0d0*tmpfac(:,8) 1624 temp(:,6)=tmpfac(:,4) 1625 temp(:,1)=tmpfac(:,5) 1626 temp(:,2)=tmpfac(:,9) 1627 end if 1628 ! Decompaction temp -> tmpfac 1629 else 1630 tmpfac(:,1:10)=0.0d0 1631 if(idir==1)then 1632 tmpfac(:,1)=3.0d0*temp(:,1) 1633 tmpfac(:,2:4)=temp(:,2:4) 1634 tmpfac(:,5:6)=2.0d0*temp(:,5:6) 1635 else if(idir==2)then 1636 tmpfac(:,2)=2.0d0*temp(:,6) 1637 tmpfac(:,9)=2.0d0*temp(:,4) 1638 tmpfac(:,4)=temp(:,5) 1639 tmpfac(:,6)=temp(:,1) 1640 tmpfac(:,8)=temp(:,3) 1641 tmpfac(:,7)=3.0d0*temp(:,2) 1642 else if(idir==3)then 1643 tmpfac(:,10)=3.0d0*temp(:,3) 1644 tmpfac(:,3)=2.0d0*temp(:,5) 1645 tmpfac(:,8)=2.0d0*temp(:,4) 1646 tmpfac(:,4)=temp(:,6) 1647 tmpfac(:,5)=temp(:,1) 1648 tmpfac(:,9)=temp(:,2) 1649 end if 1650 end if 1651 1652 end if 1653 1654 end subroutine ddkten
ABINIT/m_nonlop_pl [ Modules ]
NAME
nonlop_pl
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, GZ, MT, FF, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_nonlop_pl 22 23 use defs_basis 24 use m_errors 25 use m_abicore 26 use m_xmpi 27 use m_contistr01 28 use m_contistr03 29 use m_contistr12 30 use m_contstr21 31 use m_contstr23 32 use m_contstr25 33 use m_contstr25a 34 use m_contstr26 35 use m_metstr 36 use m_opernl 37 38 use defs_abitypes, only : MPI_type 39 use m_geometry, only : strconv 40 use m_kg, only : ph1d3d 41 use m_contract, only : cont22cso, cont22so, cont24, cont33cso, cont33so, cont35, cont22, cont3, cont13, & 42 metcon, metcon_so, metric_so 43 implicit none 44 45 private
ABINIT/nonlop_pl [ Functions ]
NAME
nonlop_pl
FUNCTION
* Compute application of a nonlocal operator Vnl in order to get: - contracted elements (energy, forces, stresses, ...), if signs=1 - a function in reciprocal space (|out> = Vnl|in>), if signs=2 Operator Vnl, as the following form: $Vnl=sum_{R,lmn,l''m''n''} {|P_{Rlmn}> Enl^{R}_{lmn,l''m''n''} <P_{Rl''m''n''}|}$ Operator Vnl is -- in the typical case -- the nonlocal potential. - With norm-conserving pseudopots, $Enl^{R}_{lmn,l''m''n''}$ is the Kleinmann-Bylander energy $Ekb^{R}_{ln}$. - The |P_{Rlmn}> are the projector functions. * This routine uses Legendre polynomials Pl to express Vnl.
INPUTS
choice: chooses possible output: choice=1 => a non-local energy contribution =2 => a gradient with respect to atomic position(s) =3 => a gradient with respect to strain(s) =23=> a gradient with respect to atm. pos. and strain(s) =4 => a gradient and 2nd derivative with respect to atomic pos. =5 => a gradient with respect to k wavevector =6 => 2nd derivatives with respect to strain dimekb1,dimekb2=dimensions of ekb (see ekb) dimffnlin=second dimension of ffnlin (1+number of derivatives) dimffnlout=second dimension of ffnlout (1+number of derivatives) ekb(dimekb1,dimekb2,nspinortot**2)= (Real) Kleinman-Bylander energies (hartree) dimekb1=lmnmax - dimekp2=ntypat ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used for the application of the nonlocal operator to the |in> vector ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used for the application of the nonlocal operator to the |out> vector gmet(3,3)=metric tensor for G vecs (in bohr**-2) gprimd(3,3)=dimensional reciprocal space primitive translations (bohr^-1) idir=direction of the - atom to be moved in the case (choice=2,signs=2), - k point direction in the case (choice=5,signs=2) - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1) indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln istwf_k=option parameter that describes the storage of wfs kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector kptin(3)=k point in terms of recip. translations, for the |in> vector kptout(3)=k point in terms of recip. translations, for the |out> vector lmnmax=max. number of (l,m,n) components over all types of atoms matblk=dimension of the arrays ph3din and ph3dout mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization natom=number of atoms in cell nattyp(ntypat)=number of atoms of each type ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpgin,nkpgout=second sizes of arrays kpgin/kpgout nloalg(3)=governs the choice of the algorithm for nonlocal operator nnlout=dimension of enlout: choice=1=>nnlout=1 choice=2=>nnlout=3*natom choice=3=>nnlout=6 choice=4=>nnlout=6*natom choice=5=>nnlout=1 choice=6=>nnlout=6*(3*natom+6) choice=23=>nnlout=6+3*natom npwin=number of planewaves for given k point, for the |in> vector npwout=number of planewaves for given k point, for the |out> vector nspinor=number of spinorial components of the wavefunctions on current proc nspinortot=total number of spinorial components of the wavefunctions ntypat=number of types of atoms in cell only_SO=flag to calculate only the SO part in nonlop phkxredin(2,natom)=phase factors exp(2 pi kptin.xred) phkxredout(2,natom)=phase factors exp(2 pi kptout.xred) ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in) ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out) --- pspso removed in beautification because it was unused --- pspso(ntypat)=spin-orbit characteristic for each atom type ------------------------------------------------------------- signs= if 1, get contracted elements (energy, forces, stress, ...) if 2, applies the non-local operator to a function in reciprocal space ucvol=unit cell volume (bohr^3) vectin(2,nspinor*npwin)=input cmplx wavefunction coefficients <G|Cnk>
OUTPUT
==== if (signs==1) ==== enlout(nnlout)= contribution of this state to the nl part of the following properties: if choice=1 : enlout(1) -> the energy if choice=2 : enlout(1:3*natom) -> the forces if choice=3 : enlout(1:6) -> the stresses if choice=23: enlout(1:6+3*natom) -> the forces and the stresses if choice=4 : enlout(1:6*natom) -> the frozen wf part of dynam. matrix if choice=6 : enlout(1:6*(3*natom+6)) -> the frozen wf part of elastic tensor ==== if (signs==2) ==== vectout(2,nspinor*npwout)= result of the aplication of the nl operator or one of its derivative to the input vect.: if choice=1 : Vnl |vectin> if choice=2 : dVnl/d(xred(idir,iatom) |vectin> (xred=reduced atm. pos.) if choice=3 : dVnl/d(strain(idir)) |vectin> (symmetric strain =>idir=1...6) if choice=5 : dVnl/dk(idir) |vectin> (k wavevector)
NOTES
In the case signs=1, the array vectout is not used, nor modified so that the same array as vectin can be used as a dummy argument; the same is true for the pairs npwin-npwout, ffnlin-ffnlout, kgin-kgout, ph3din-ph3dout, phkredin-phkxredout). Calculation includes contributions to grads of Etot wrt coord and wrt strains for l=0,1,2,3.
WARNINGS
- Warning 1: This routine is in a transient state, during the time of the implementation of the spin-orbit coupling... In particular, the OMP parallelisation is still missing, but it matters here only when nspinor==2. - Warning 2: the order of atoms is governed by atindx
SOURCE
170 subroutine nonlop_pl(choice,dimekb1,dimekb2,dimffnlin,dimffnlout,ekb,enlout,& 171 & ffnlin,ffnlout,gmet,gprimd,idir,indlmn,istwf_k,kgin,kgout,kpgin,kpgout,& 172 & kptin,kptout,lmnmax,matblk,mgfft,mpi_enreg,mpsang,mpssoang,& 173 & natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,npwin,npwout,nspinor,nspinortot,& 174 & ntypat,only_SO,phkxredin,phkxredout,ph1d,ph3din,ph3dout,signs,& 175 & ucvol,vectin,vectout) 176 177 !Arguments ------------------------------------ 178 !This type is defined in defs_mpi 179 !The (inout) classification below is misleading; mpi_enreg is temporarily 180 ! changed but reset to its initial condition before exiting. 181 !scalars 182 integer,intent(in) :: choice,dimekb1,dimekb2,dimffnlin,dimffnlout,idir,istwf_k 183 integer,intent(in) :: lmnmax,matblk,mgfft,mpsang,mpssoang,natom,nkpgin,nkpgout 184 integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,only_SO,signs 185 real(dp),intent(in) :: ucvol 186 type(MPI_type),intent(in) :: mpi_enreg 187 !arrays 188 integer,intent(in) :: indlmn(6,lmnmax,ntypat),kgin(3,npwin),kgout(3,npwout) 189 integer,intent(in) :: nattyp(ntypat),ngfft(18),nloalg(3) !,pspso(ntypat) UNUSED 190 real(dp),intent(in) :: ekb(dimekb1,dimekb2,nspinortot**2) 191 real(dp),intent(in) :: ffnlin(npwin,dimffnlin,lmnmax,ntypat) 192 real(dp),intent(in) :: ffnlout(npwout,dimffnlout,lmnmax,ntypat),gmet(3,3) 193 real(dp),intent(in) :: gprimd(3,3),kpgin(npwin,nkpgin),kpgout(npwout,nkpgout) 194 !real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom) !vz_d 195 real(dp),intent(in) :: kptin(3),kptout(3) !vz_d 196 real(dp),intent(in) :: ph1d(2,*) !vz_d 197 real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom) 198 real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk) 199 real(dp),intent(inout) :: vectin(:,:) 200 real(dp),intent(out) :: enlout(:) !vz_i 201 real(dp),intent(inout) :: vectout(:,:) !vz_i 202 203 !Local variables------------------------------- 204 !mlang is the maximum number of different angular momenta 205 !(mlang=4 means s,p,d,f) 206 ! Note : in a future version, one should adjust mlang to mpsang. 207 !mlang2 is the maximum number of unique tensor components for a tensor 208 !of rank (mlang-1) with index range 1-3 209 !mlang3 is the maximum number of unique tensor components summed over 210 !all tensors of rank 0 through mlang-1. 211 !mlang4 is the total number of additional unique tensor components 212 !related to strain gradients, ranks 2 through mlang+1. 213 !mlang6 is the total number of additional unique tensor components 214 !related to strain 2nd derivaives, ranks 4 through mlang+3. 215 !mlang1 is the total number of certain additional unique tensor components 216 !related to internal strain, ranks 1 through mlang 217 !mlang5 is the total number of other additional unique tensor components 218 !related to internal strain, ranks 1 through mlang 219 !scalars 220 integer,parameter :: mlang=4 221 ! MG: I tried to use parameters instead of saved variables but [tutorespfn][trf2_1] gets stuck on milou_g95_snofbfpe 222 integer,save :: mlang1=((mlang+1)*(mlang+2)*(mlang+3))/6-1 223 !integer,save :: mlang2=(mlang*(mlang+1))/2 ! Unused 224 integer,save :: mlang3=(mlang*(mlang+1)*(mlang+2))/6 225 integer,save :: mlang4=((mlang+2)*(mlang+3)*(mlang+4))/6-4 226 integer,save :: mlang5=((mlang+3)*(mlang+4)*(mlang+5))/6-10 227 integer,save :: mlang6=((mlang+4)*(mlang+5)*(mlang+6))/6-20 228 integer :: compact,ia,ia1,ia2,ia3,ia4,ia5,ierr,iest,ig,ii,ilang,ilang2,ilmn 229 integer :: iln,iln0,indx,iproj,ipsang,ishift,isp,ispin,ispinor,ispinor_index,ispinp 230 integer :: istr,istr1,istr2,iterm,itypat,jj,jjk,jjs,jjs1,jjs2,jjs3,jjs4,jjstr,jspin 231 integer :: mincat,mproj,mu,mumax,n1,n2,n3,ndgxdt,ndgxdtfac,nincat,nlang 232 integer :: nproj,nspinso,rank 233 integer :: sign,spaceComm, isft 234 real(dp) :: e2nl,e2nldd,enlk 235 character(len=500) :: msg 236 !arrays 237 integer,allocatable :: indlmn_s(:,:,:),jproj(:) 238 real(dp) :: amet(2,3,3,2,2),amet_lo(3,3),e2nl_tmp(6),eisnl(3),rank2(6) 239 real(dp) :: rank2c(2,6),strsnl(6),strsnl_out(6),strsso(6,3),strssoc(6),trace(2)!,tsec(2) 240 real(dp),allocatable :: d2gxdis(:,:,:,:,:),d2gxdis_s(:,:,:,:) 241 real(dp),allocatable :: d2gxds2(:,:,:,:,:),d2gxds2_s(:,:,:,:) 242 real(dp),allocatable :: dgxdis(:,:,:,:,:),dgxdis_s(:,:,:,:),dgxds(:,:,:,:,:) 243 real(dp),allocatable :: dgxds_s(:,:,:,:),dgxdsfac(:,:,:,:,:) 244 real(dp),allocatable :: dgxdt(:,:,:,:,:,:),dgxdt_s(:,:,:,:,:) 245 real(dp),allocatable :: dgxdtfac(:,:,:,:,:),ekb_s(:,:),gxa(:,:,:,:,:) 246 real(dp),allocatable :: gxa_s(:,:,:,:),gxafac(:,:,:,:),pauli(:,:,:,:) 247 real(dp),allocatable :: temp(:,:),tmpfac(:,:),vectin_s(:,:),vectout_s(:,:) 248 real(dp),allocatable :: wt(:,:) 249 250 ! ********************************************************************** 251 252 ABI_UNUSED(mgfft) 253 254 !Test: spin orbit not allowed for choice=5,6 255 if (nspinortot==2 .and. choice==6 ) then 256 ABI_BUG('For nspinortot=2, choice=6 is not yet allowed.') 257 end if 258 259 if ((choice<1 .or. choice>6) .and. choice/=23 ) then 260 write(msg,'(a,i0)')'Does not presently support this choice=',choice 261 ABI_BUG(msg) 262 end if 263 264 !Test: choice 51 and 52 only allowed with nonlop_ylm 265 !JWZ, 01-Sep-08 266 if (choice==51 .or. choice==52) then 267 ABI_BUG('choice 51 or 52 is not yet allowed.') 268 end if 269 270 !Define dimension of work arrays. 271 mincat=min(NLO_MINCAT,maxval(nattyp)) 272 mproj=maxval(indlmn(3,:,:)) 273 ABI_MALLOC(temp,(2,mlang4)) 274 ABI_MALLOC(tmpfac,(2,mlang4)) 275 ABI_MALLOC(wt,(mlang,mproj)) 276 ABI_MALLOC(jproj,(mlang)) 277 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 278 ABI_MALLOC(ekb_s,(mlang,mproj)) 279 ABI_MALLOC(indlmn_s,(6,lmnmax,ntypat)) 280 281 !Eventually compute the spin-orbit metric tensor: 282 if (mpssoang>mpsang) then 283 ABI_MALLOC(pauli,(2,2,2,3)) 284 call metric_so(amet,gprimd,pauli) 285 end if 286 287 !Allocate array gxa (contains projected scalars). 288 ABI_MALLOC(gxa,(2,mlang3,mincat,mproj,nspinortot)) 289 if(nspinor==2) then 290 ABI_MALLOC(gxa_s,(2,mlang3,mincat,mproj)) 291 else 292 ABI_MALLOC(gxa_s,(0,0,0,0)) 293 end if 294 295 ABI_MALLOC(gxafac,(2,mlang3,mincat,mproj)) 296 gxa(:,:,:,:,:)=zero 297 298 !If choice==2 : first-order atomic displacements 299 !If signs==2, only one direction is considered 300 !If signs==1, the three directions are considered 301 !If choice==4 and signs==1 : second-order atomic displacements, 302 !the nine components are considered 303 !If choice==5 and signs==2 : ddk 304 !component 1 -> from ffnl(:,2,...) 305 !component 2 -> from ffnl(:,1,...) (too much space is booked for this 306 !case, since the number of angular momenta is smaller than mlang3, but it is easier) 307 ndgxdt=0 308 if(signs==2 .and. choice==2) ndgxdt=1 309 if(signs==1 .and. (choice==2.or.choice==23)) ndgxdt=3 310 if(choice==4) ndgxdt=9 311 if(choice==5) ndgxdt=2 312 !Allocate dgxdt (contains derivatives of gxa with respect to atomic displacements or ddk). 313 ABI_MALLOC(dgxdt,(2,ndgxdt,mlang3,mincat,mproj,nspinortot)) 314 dgxdt(:,:,:,:,:,:)=zero 315 if(nspinor==2)then 316 ABI_MALLOC(dgxdt_s,(2,ndgxdt,mlang3,mincat,mproj)) 317 dgxdt_s(:,:,:,:,:)=zero 318 else 319 ABI_MALLOC(dgxdt_s,(0,0,0,0,0)) 320 end if 321 ndgxdtfac=0 322 if(signs==2 .and. choice==2) ndgxdtfac=1 323 if(choice==4) ndgxdtfac=3 324 if(choice==5) ndgxdtfac=2 325 ABI_MALLOC(dgxdtfac,(2,ndgxdtfac,mlang3,mincat,mproj)) 326 327 !Allocate dgxds (contains derivatives of gxa with respect to strains). 328 ABI_MALLOC(dgxds,(2,mlang4,mincat,mproj,nspinor)) 329 dgxds(:,:,:,:,:)=zero 330 ABI_MALLOC(dgxdsfac,(2,mlang4,mincat,mproj,nspinor)) 331 if(choice==6) then 332 ABI_MALLOC(dgxdis,(2,mlang1,mincat,mproj,nspinor)) 333 ABI_MALLOC(d2gxdis,(2,mlang5,mincat,mproj,nspinor)) 334 ABI_MALLOC(d2gxds2,(2,mlang6,mincat,mproj,nspinor)) 335 else 336 ABI_MALLOC(dgxdis ,(0,0,0,0,0)) 337 ABI_MALLOC(d2gxdis,(0,0,0,0,0)) 338 ABI_MALLOC(d2gxds2,(0,0,0,0,0)) 339 end if 340 ABI_MALLOC(dgxds_s ,(0,0,0,0)) 341 ABI_MALLOC(dgxdis_s ,(0,0,0,0)) 342 ABI_MALLOC(d2gxdis_s,(0,0,0,0)) 343 ABI_MALLOC(d2gxds2_s,(0,0,0,0)) 344 if(nspinor==2)then 345 ABI_FREE(dgxds_s) 346 ABI_MALLOC(dgxds_s,(2,mlang4,mincat,mproj)) 347 dgxds_s(:,:,:,:)=zero 348 if(choice==6) then 349 ABI_FREE(dgxdis_s) 350 ABI_FREE(d2gxdis_s) 351 ABI_FREE(d2gxds2_s) 352 ABI_MALLOC(dgxdis_s,(2,mlang1,mincat,mproj)) 353 ABI_MALLOC(d2gxdis_s,(2,mlang5,mincat,mproj)) 354 ABI_MALLOC(d2gxds2_s,(2,mlang6,mincat,mproj)) 355 else 356 end if 357 end if 358 359 !Zero out some arrays 360 if(choice==2 .or. choice==4 .or. choice==5 .or. choice==6 .or. choice==23) enlout(:)=0.0d0 361 362 if(signs==2) vectout(:,:)=zero 363 364 !if(choice==3.or.choice==23) then 365 ! strsnl(:)=zero 366 ! if(mpssoang>mpsang) strsso(:,:)=zero 367 !end if 368 enlk=zero 369 strsso = zero 370 strsnl = zero 371 372 !In the case vectin is a spinor, split its second part. 373 !Also, eventually take into account the storage format of the wavefunction 374 !(the original vector will be restored before leaving the routine, 375 !except for the vectin(2,1) component with istwf_k==2, 376 !that should vanish) 377 !In sequential, treat the second spinor part first 378 if (nspinor==2)then 379 ABI_MALLOC(vectin_s,(2,npwin)) 380 ABI_MALLOC(vectout_s,(2,npwout)) 381 382 isft = npwin;if (mpi_enreg%nproc_spinor>1) isft=0 383 384 ! Initialize it 385 !$OMP PARALLEL DO 386 do ig=1,npwin 387 vectin_s(1,ig)=vectin(1,ig+isft) 388 vectin_s(2,ig)=vectin(2,ig+isft) 389 end do 390 391 ! Take into account the storage 392 if(istwf_k/=1)then 393 call scalewf_nonlop(istwf_k,mpi_enreg,npwin,1,vectin_s) 394 end if 395 end if ! nspinortot==2 396 397 !Treat the first spinor part now 398 if(istwf_k/=1) then 399 call scalewf_nonlop(istwf_k,mpi_enreg,npwin,1,vectin) 400 end if 401 402 403 !Big loop on atom types. 404 ia1=1 405 do itypat=1,ntypat 406 407 ! Get atom loop indices for different types: 408 ia2=ia1+nattyp(itypat)-1 409 410 ! Cut the sum on different atoms in blocks, to allow memory saving. 411 ! Inner summations on atoms will be done from ia3 to ia4. 412 ! Note that the maximum range from ia3 to ia4 is mincat (maximum 413 ! increment of atoms). 414 do ia3=ia1,ia2,mincat 415 ia4=min(ia2,ia3+mincat-1) 416 ! Give the increment of number of atoms in this subset. 417 nincat=ia4-ia3+1 418 419 ! Prepare the phase factors for the atoms between ia3 and ia4 : 420 ! For nloalg(2)<=0, they were not prepared previously,it is needed to 421 ! compute them again. 422 if(nloalg(2)<=0)then 423 ! For nloalg(2)==0, it is needed to compute the phase factors. 424 if(mincat>matblk)then 425 write(msg,'(a,a,a,i4,a,i4,a)')& 426 'With nloc_mem<=0, mincat must be less than matblk.',ch10,& 427 'Their value is ',mincat,' and ',matblk,'.' 428 ABI_BUG(msg) 429 end if 430 call ph1d3d(ia3,ia4,kgin,matblk,natom,npwin,n1,n2,n3,phkxredin,ph1d,ph3din) 431 end if 432 433 ! Here begins the different treatment for the scalar-relativistic 434 ! part(s) and the spin-orbit part. 435 ! Loop on ispinor : 1 for scalar-Relativistic, 2 for spin-orbit 436 nspinso=1;if (mpssoang>mpsang) nspinso=2 437 438 ! Change nspinso if collinear run or if nspinor == 2 and SOC is not wanted. 439 ! TODO: The last check requires pspso 440 if (nspinortot == 1) nspinso = 1 441 442 do ispinor=1,nspinso 443 ispinor_index=ispinor 444 if (mpi_enreg%paral_spinor==1) ispinor_index=mpi_enreg%me_spinor+1 445 446 ! 447 ! mjv 25 6 2008: if only_SO == 1 or 2 skip the scalar relativistic terms 448 ! only output the spin orbit ones 449 ! 450 if (ispinor==1 .and. only_SO>0) cycle 451 452 ! Select scalar-relativistic or spin-orbit KB-energies: 453 ekb_s(:,:)=zero ; wt(:,:)=zero 454 ! Loop over (l,n) values (loop over l,m,n and test on l,n) 455 iln0=0 ; jproj(:)=0 ; nlang=0 456 indlmn_s(:,:,itypat)=0 457 do ilmn=1,lmnmax 458 if(ispinor/=indlmn(6,ilmn,itypat))cycle 459 indlmn_s(:,ilmn,itypat)=indlmn(:,ilmn,itypat) 460 iln =indlmn(5,ilmn,itypat) 461 if (iln>iln0) then 462 iln0=iln 463 ipsang=indlmn(1,ilmn,itypat)+1 464 ! DEBUG 465 ! write(std_out,*)' nonlop : ipsang,ilmn,itypat,ispinor=',ipsang,ilmn,itypat,ispinor 466 ! ENDDEBUG 467 iproj=indlmn(3,ilmn,itypat) 468 ! This shift is not needed anymore 469 ! if (ispinor==2) ipsang=indlmn(1,ilmn,itypat)-mpsang+2 470 ekb_s(ipsang,iproj)=ekb(iln,itypat,ispinor) 471 wt(ipsang,iproj)=4.d0*pi/ucvol*dble(2*ipsang-1)*ekb_s(ipsang,iproj) 472 ! 473 ! mjv 27 6 2008: if only_SO == 2 remove the factor of l in the operator 474 ! 475 if (only_SO == 2) then 476 wt(ipsang,iproj)=4.d0*pi/ucvol*ekb_s(ipsang,iproj) 477 end if 478 jproj(ipsang)=max(jproj(ipsang),iproj) 479 if(iproj>0)nlang=max(nlang,ipsang) 480 end if 481 end do ! ilmn 482 483 484 ! If nlang is not 0, then some non-local part is to be applied for that type of atom. 485 if (nlang/=0) then 486 ! Operate with the non-local potential on the wavefunction, in order 487 ! to get projected quantities. Call different routines opernl2, 488 ! opernl3, opernl4x which corresponds to different writings of the 489 ! same numerical operations. There is still optimisation left for 490 ! the case istwf_k/=1 (up to a factor 2 in speed). 491 ! call timab(74+choice,1,tsec) 492 sign=1 493 gxa(:,:,:,:,:)=zero 494 dgxdt(:,:,:,:,:,:)=zero 495 dgxds(:,:,:,:,:)=zero 496 497 ! Only the first spinorial component of vectin is taken into account first 498 ispin=1;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index 499 if(nloalg(1)==2) then 500 call opernl2(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 501 & ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 502 & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 503 & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,& 504 & ntypat,ph3din,sign,vectin) 505 else if(nloalg(1)==3) then 506 call opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 507 & ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 508 & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 509 & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,& 510 & ntypat,ph3din,sign,vectin) 511 else if(nloalg(1)==4) then 512 call opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 513 & ffnlin,gmet,gxa,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 514 & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 515 & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,& 516 & ntypat,ph3din,vectin) 517 end if 518 ! This duplication of the opernl calls is needed to avoid copying 519 ! vectin, with a detrimental effect on speed. 520 if (nspinor==2)then 521 ispin=2 522 if(nloalg(1)==2) then 523 call opernl2(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,& 524 & ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 525 & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 526 & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,& 527 & ntypat,ph3din,sign,vectin_s) 528 else if(nloalg(1)==3) then 529 call opernl3(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,& 530 & ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 531 & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 532 & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,& 533 & ntypat,ph3din,sign,vectin_s) 534 else if(nloalg(1)==4) then 535 call opernl4a(choice,dgxdis_s,dgxds_s,d2gxdis_s,d2gxds2_s,dgxdt_s,& 536 & ffnlin,gmet,gxa_s,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 537 & jproj,kgin,kpgin,kptin,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 538 & mlang5,mlang6,mproj,ndgxdt,dimffnlin,nincat,nkpgin,nlang,nloalg,npwin,& 539 & ntypat,ph3din,vectin_s) 540 end if 541 dgxds(:,:,:,:,ispin)=dgxds_s(:,:,:,:) 542 dgxdt(:,:,:,:,:,ispin)=dgxdt_s(:,:,:,:,:) 543 gxa(:,:,:,:,ispin)=gxa_s(:,:,:,:) 544 end if 545 546 ! Parallelism stuff 547 spaceComm=mpi_enreg%comm_fft 548 call xmpi_sum(dgxds,spaceComm,ierr) 549 if (mpi_enreg%paral_spinor==1) then 550 spaceComm=mpi_enreg%comm_spinorfft 551 jspin=3-ispin 552 gxa(:,:,:,:,ispin)=gxa(:,:,:,:,1) 553 gxa(:,:,:,:,jspin)=zero 554 if ( ndgxdt>0)then 555 dgxdt(:,:,:,:,:,ispin)=dgxdt(:,:,:,:,:,1) 556 dgxdt(:,:,:,:,:,jspin)=zero 557 end if 558 end if 559 560 call xmpi_sum(gxa,spaceComm,ierr) 561 if(ndgxdt>0) then 562 call xmpi_sum(dgxdt,spaceComm,ierr) 563 end if 564 565 ! XG030513 : MPIWF, at this place, one should perform the reduction 566 ! and spread of data of gxa, dgxdt and dgxds 567 568 ! Loop over spins: 569 do isp=1,nspinor 570 ispin=isp;if (mpi_enreg%paral_spinor==1) ispin=ispinor_index 571 572 ! Perform contractions for the various tensors (d)gx?, producing the 573 ! contracted tensors (d)gx?fac to be passed back to opernl: 574 do ia=1,nincat 575 do ilang=1,nlang 576 nproj=jproj(ilang) 577 if(nproj/=0) then 578 ilang2=(ilang*(ilang+1))/2 579 do iproj=1,nproj 580 581 ! The rank of the tensor gxa equals l: 582 rank=ilang-1 583 ! jjs gives the starting address of the relevant components 584 jjs=1+((ilang-1)*ilang*(ilang+1))/6 585 if (ilang>4) then 586 write(msg,'(a,i0)')' ilang must fall in range [1..4] but value is ',ilang 587 ABI_BUG(msg) 588 end if 589 590 ! Metric & spinorial contraction from gxa to gxafac. The treatment 591 ! is different for the scalar-relativistic and spin-orbit parts. 592 if(ispinor==1) then 593 ! ------ Scalar-Relativistic ------ 594 temp(:,1:((rank+1)*(rank+2))/2)= & 595 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 596 call metcon(rank,gmet,temp,tmpfac) 597 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= & 598 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 599 else 600 ! ------ Spin-orbit ------ 601 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero 602 ! Contraction over spins: 603 do ispinp=1,nspinortot 604 ! => Imaginary part (multiplying by i, then by the Im of amet): 605 temp(1,1:((rank+1)*(rank+2))/2)= & 606 & -gxa(2,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp) 607 temp(2,1:((rank+1)*(rank+2))/2)= & 608 & gxa(1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp) 609 amet_lo(:,:)=amet(2,:,:,ispin,ispinp) 610 call metcon_so(rank,gmet,amet_lo,temp,tmpfac) 611 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= & 612 & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ & 613 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 614 ! => Real part: 615 temp(:,1:((rank+1)*(rank+2))/2)= & 616 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp) 617 amet_lo(:,:)=amet(1,:,:,ispin,ispinp) 618 call metcon_so(rank,gmet,amet_lo,temp,tmpfac) 619 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= & 620 & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ & 621 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 622 end do 623 end if 624 625 ! Eventual tensorial compaction/decompaction for ddk 626 ! perturbation: 627 if(choice==5 .and. ilang>=2) then 628 jjk=1+((ilang-2)*(ilang-1)*ilang)/6 629 compact=-1 630 temp(:,1:(rank*(rank+1))/2)= & 631 & dgxdt(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj,ispin) 632 call ddkten(compact,idir,rank,temp,tmpfac) 633 dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)= & 634 & dgxdt(:,1,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)& 635 & +tmpfac(:,1:((rank+1)*(rank+2))/2) 636 compact=1 637 tmpfac(:,1:((rank+1)*(rank+2))/2)= & 638 & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj) 639 call ddkten(compact,idir,rank,temp,tmpfac) 640 dgxdtfac(:,2,jjk:jjk-1+(rank*(rank+1))/2,ia,iproj)= & 641 & temp(:,1:(rank*(rank+1))/2) 642 end if 643 644 ! Section for strain perturbation 645 ! no spin-orbit yet 646 647 if(choice==3 .and. signs==2) then 648 istr=idir 649 if(ispinor==1) then 650 ! ------ Scalar-Relativistic ------ 651 ! jjstr is the starting address for dgxds and dgxdsfac 652 jjstr=-3+((ilang+1)*(ilang+2)*(ilang+3))/6 653 ! diagonal enlk contribution 654 ! note sign change (12/05/02) 655 if(istr>3) then 656 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero 657 else 658 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=& 659 & -gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj) 660 end if 661 dgxdsfac(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)=zero 662 iterm=1 663 temp(:,1:((rank+1)*(rank+2))/2)= & 664 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 665 call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac) 666 dgxdsfac(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp)= & 667 & wt(ilang,iproj)*tmpfac(:,1:((rank+3)*(rank+4))/2) 668 iterm=2 669 temp(:,1:((rank+1)*(rank+2))/2)= & 670 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 671 call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac) 672 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= & 673 & +gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ & 674 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 675 iterm=3 676 temp(:,1:((rank+3)*(rank+4))/2)= & 677 dgxds(:,jjstr:jjstr-1+((rank+3)*(rank+4))/2,ia,iproj,isp) 678 call metstr(istr,rank,iterm,gmet,gprimd,temp,tmpfac) 679 gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)= & 680 & gxafac(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+ & 681 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 682 end if 683 ! end section for strain perturbation 684 end if 685 686 ! Eventual metric & spinorial contraction from dgxdt to dgxdtfac. 687 ! either for the dynamical matrix, or for the application of the 688 ! gradient of the operator. The treatment is different for the 689 ! scalar-relativistic and spin-orbit parts. 690 if ((choice==2.and.signs==2).or. & 691 & (choice==5.and.signs==2).or. & 692 & (choice==4)) then 693 mumax=ndgxdtfac;if (choice==5) mumax=1 694 do mu=1,mumax 695 if(ispinor==1) then 696 ! ------ Scalar-Relativistic ------ 697 temp(:,1:((rank+1)*(rank+2))/2)= & 698 & dgxdt(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 699 call metcon(rank,gmet,temp,tmpfac) 700 dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=& 701 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 702 else 703 ! ------ Spin-orbit ------ 704 dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=zero 705 ! Contraction over spins: 706 do ispinp=1,nspinortot 707 ! => Imaginary part (multiplying by i, then by the Im of amet): 708 temp(1,1:((rank+1)*(rank+2))/2)= & 709 & -dgxdt(2,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp) 710 temp(2,1:((rank+1)*(rank+2))/2)= & 711 & dgxdt(1,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp) 712 amet_lo(:,:)=amet(2,:,:,ispin,ispinp) 713 call metcon_so(rank,gmet,amet_lo,temp,tmpfac) 714 dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=& 715 & dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+& 716 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 717 ! => Real part: 718 temp(:,1:((rank+1)*(rank+2))/2)= & 719 & dgxdt(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispinp) 720 amet_lo(:,:)=amet(1,:,:,ispin,ispinp) 721 call metcon_so(rank,gmet,amet_lo,temp,tmpfac) 722 dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)=& 723 & dgxdtfac(:,mu,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj)+& 724 & wt(ilang,iproj)*tmpfac(:,1:((rank+1)*(rank+2))/2) 725 end do 726 end if 727 end do 728 end if 729 730 731 ! ---- Accumulate the nonlocal energy. 732 do ii=1,ilang2 733 jj=ii-1+jjs 734 enlk=enlk+(gxafac(1,jj,ia,iproj)*gxa(1,jj,ia,iproj,ispin)& 735 & +gxafac(2,jj,ia,iproj)*gxa(2,jj,ia,iproj,ispin) ) 736 end do 737 738 ! ---- Accumulate force contributions if requested. 739 ! Note that the contraction of gxa with dgxdt can be done like 740 ! a cartesian dot product now because the symmetrically-related 741 ! terms are accounted for with weights in gxa. 742 if ((choice==2.or.choice==23) .and. signs==1) then 743 ishift=0;if (choice==23) ishift=6 744 ia5=ia+ia3-1 745 do ii=1,ilang2 746 jj=ii-1+jjs 747 do mu=1,3 748 ! (includes also factor of 2 from "2*Re[stuff]") 749 indx=mu+3*(ia5-1)+ishift 750 enlout(indx)=enlout(indx)+two*& 751 & ( gxafac(1,jj,ia,iproj)*dgxdt(1,mu,jj,ia,iproj,ispin)& 752 & +gxafac(2,jj,ia,iproj)*dgxdt(2,mu,jj,ia,iproj,ispin)) 753 end do 754 end do 755 end if 756 757 ! ---- Accumulate stress tensor contributions if requested. 758 if ((choice==3.or.choice==23).and.signs==1) then 759 ! 1- Compute contractions involving gxa and dgxds: 760 ! --- Same formula for Scalar-relativistic and Spin-orbit --- 761 if (ilang==1) then 762 do ii=1,6 763 rank2(ii)=2.d0*& 764 & (gxafac(1,1,ia,iproj)*dgxds(1,ii,ia,iproj,isp)+ & 765 & gxafac(2,1,ia,iproj)*dgxds(2,ii,ia,iproj,isp) ) 766 end do 767 else if (ilang==2) then 768 call cont13(gxafac(:,jjs:jjs+2,ia,iproj),& 769 & dgxds(:, 7:16,ia,iproj,isp),rank2) 770 else if (ilang==3) then 771 call cont24(gxafac(:,jjs:jjs+5,ia,iproj),& 772 & dgxds(:,17:31,ia,iproj,isp),rank2) 773 else if (ilang==4) then 774 call cont35(gxafac(:,jjs:jjs+9,ia,iproj),& 775 & dgxds(:,32:52,ia,iproj,isp),rank2) 776 end if 777 ! In all cases add rank2 term into stress tensor 778 strsnl(:)=strsnl(:)-rank2(:) 779 ! 2- Compute contractions involving gxa and gxa: 780 if(ispinor==1) then 781 ! 2a ------ Scalar-Relativistic ------ 782 if (ilang==2) then 783 strsnl(1)=strsnl(1)-wt(ilang,iproj)*2.d0*& 784 & (gxa(1,jjs ,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispin)+& 785 & gxa(2,jjs ,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispin)) 786 strsnl(2)=strsnl(2)-wt(ilang,iproj)*2.d0*& 787 & (gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispin)+& 788 & gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispin)) 789 strsnl(3)=strsnl(3)-wt(ilang,iproj)*2.d0*& 790 & (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+2,ia,iproj,ispin)+& 791 & gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+2,ia,iproj,ispin)) 792 strsnl(4)=strsnl(4)-wt(ilang,iproj)*2.d0*& 793 & (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispin)+& 794 & gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispin)) 795 strsnl(5)=strsnl(5)-wt(ilang,iproj)*2.d0*& 796 & (gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispin)+& 797 & gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispin)) 798 strsnl(6)=strsnl(6)-wt(ilang,iproj)*2.d0*& 799 & (gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispin)+& 800 & gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispin)) 801 else if (ilang==3) then 802 call trace2(gxa(:,jjs:jjs+5,ia,iproj,ispin),gmet,trace) 803 call cont22(gxa(:,jjs:jjs+5,ia,iproj,ispin),gmet,rank2) 804 do ii=1,6 805 strsnl(ii)=strsnl(ii)+wt(ilang,iproj)*& 806 & (2.d0*(trace(1)*gxa(1,jjs+ii-1,ia,iproj,ispin)+& 807 & trace(2)*gxa(2,jjs+ii-1,ia,iproj,ispin))-3.d0*rank2(ii)) 808 end do 809 else if (ilang==4) then 810 call cont3(gxa(:,jjs:jjs+9,ia,iproj,ispin),gmet,rank2) 811 strsnl(:)=strsnl(:)-wt(ilang,iproj)*rank2(:) 812 end if 813 else 814 ! 2b ------ Spin-orbit ------ 815 do ispinp=1,nspinortot 816 if (ilang==3) then 817 call cont22so(gxa(:,jjs:jjs+5,ia,iproj,ispin),& 818 & gxa(:,jjs:jjs+5,ia,iproj,ispinp),& 819 & amet(:,:,:,ispin,ispinp),rank2) 820 strsnl(:)=strsnl(:)-wt(ilang,iproj)*3.d0*rank2(:) 821 else if (ilang==4) then 822 call cont33so(gxa(:,jjs:jjs+9,ia,iproj,ispin),& 823 & gxa(:,jjs:jjs+9,ia,iproj,ispinp),& 824 & gmet,amet(:,:,:,ispin,ispinp),rank2) 825 strsnl(:)=strsnl(:)-wt(ilang,iproj)*rank2(:) 826 end if 827 end do 828 end if 829 ! 3- Compute contractions involving gxa and gxa due to 830 ! gradients of antisymmetric tensor (amet): 831 ! --- Only in case of Spin-orbit --- 832 if(ispinor==2) then 833 do ispinp=1,nspinortot 834 ! Be carefull: no need to compute rank2c(:,1:3) ! 835 if (ilang==2) then 836 rank2c(1,4)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispinp)& 837 & +gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispinp) 838 rank2c(2,4)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(2,jjs+1,ia,iproj,ispinp)& 839 & -gxa(2,jjs+2,ia,iproj,ispin)*gxa(1,jjs+1,ia,iproj,ispinp) 840 rank2c(1,5)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp)& 841 & +gxa(2,jjs+2,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp) 842 rank2c(2,5)=gxa(1,jjs+2,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp)& 843 & -gxa(2,jjs+2,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp) 844 rank2c(1,6)=gxa(1,jjs+1,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp)& 845 & +gxa(2,jjs+1,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp) 846 rank2c(2,6)=gxa(1,jjs+1,ia,iproj,ispin)*gxa(2,jjs ,ia,iproj,ispinp)& 847 & -gxa(2,jjs+1,ia,iproj,ispin)*gxa(1,jjs ,ia,iproj,ispinp) 848 else if (ilang==3) then 849 call cont22cso(gxa(:,jjs:jjs+5,ia,iproj,ispin),& 850 & gxa(:,jjs:jjs+5,ia,iproj,ispinp),& 851 & gmet,rank2c) 852 else if (ilang==4) then 853 call cont33cso(gxa(:,jjs:jjs+9,ia,iproj,ispin),& 854 & gxa(:,jjs:jjs+9,ia,iproj,ispinp),& 855 & gmet,rank2c) 856 end if 857 if (ilang>1) then 858 do jj=1,3 859 do ii=4,6 860 strsso(ii,jj)=strsso(ii,jj)-2.d0*wt(ilang,iproj)*& 861 & (pauli(1,ispin,ispinp,jj)*rank2c(2,ii)+& 862 & pauli(2,ispin,ispinp,jj)*rank2c(1,ii)) 863 end do 864 end do 865 end if 866 end do 867 end if 868 ! Enf if (choice==3) 869 end if 870 871 ! ---- Accumulate dynamical matrix contributions if requested. 872 if (choice==4) then 873 ia5=ia+ia3-1 874 do ii=1,ilang2 875 jj=ii-1+jjs 876 do mu=1,6 877 ! (includes also factor of 2 from "2*Re[stuff]") 878 enlout(mu+6*(ia5-1))=enlout(mu+6*(ia5-1))+two*& 879 & (gxafac(1,jj,ia,iproj)*dgxdt(1,mu+3,jj,ia,iproj,ispin)& 880 & +gxafac(2,jj,ia,iproj)*dgxdt(2,mu+3,jj,ia,iproj,ispin)) 881 end do 882 do mu=1,3 883 enlout(mu+6*(ia5-1))=enlout(mu+6*(ia5-1))+two*& 884 & (dgxdtfac(1,mu,jj,ia,iproj)*dgxdt(1,mu,jj,ia,iproj,ispin)& 885 & +dgxdtfac(2,mu,jj,ia,iproj)*dgxdt(2,mu,jj,ia,iproj,ispin)) 886 end do 887 enlout(4+6*(ia5-1))=enlout(4+6*(ia5-1)) +two*& 888 & (dgxdtfac(1,2,jj,ia,iproj)*dgxdt(1,3,jj,ia,iproj,ispin)& 889 & +dgxdtfac(2,2,jj,ia,iproj)*dgxdt(2,3,jj,ia,iproj,ispin)) 890 enlout(5+6*(ia5-1))=enlout(5+6*(ia5-1)) +two*& 891 & (dgxdtfac(1,1,jj,ia,iproj)*dgxdt(1,3,jj,ia,iproj,ispin)& 892 & +dgxdtfac(2,1,jj,ia,iproj)*dgxdt(2,3,jj,ia,iproj,ispin)) 893 enlout(6+6*(ia5-1))=enlout(6+6*(ia5-1)) +two*& 894 & (dgxdtfac(1,1,jj,ia,iproj)*dgxdt(1,2,jj,ia,iproj,ispin)& 895 & +dgxdtfac(2,1,jj,ia,iproj)*dgxdt(2,2,jj,ia,iproj,ispin)) 896 end do 897 end if 898 899 ! ---- Accumulate elastic tensor contributions if requested. 900 901 if(choice==6) then 902 ! XG 081121 : msg to the person who has introduced this CPP option (sorry, I did not have to time to locate who did this ...) 903 ! This section of ABINIT should be allowed by default to the user. I have found that on the contrary, the build 904 ! system defaults are such that this section is forbidden by default. You might restore this flag if at the same time, 905 ! you modify the build system in such a way that by default this section is included, and if the user wants, it can disable it. 906 ! #if defined USE_NLSTRAIN_LEGENDRE 907 ia5=ia+ia3-1 908 jjs1=((ilang)*(ilang+1)*(ilang+2))/6 909 jjs2=-3+((ilang+1)*(ilang+2)*(ilang+3))/6 910 jjs3=-9+((ilang+2)*(ilang+3)*(ilang+4))/6 911 jjs4=-19+((ilang+3)*(ilang+4)*(ilang+5))/6 912 913 ! Contribution for two diagonal strains (basically, the nonlocal 914 ! enlk) 915 temp(:,1:((rank+1)*(rank+2))/2)= & 916 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 917 call metcon(rank,gmet,temp,tmpfac) 918 e2nldd=zero 919 do ii=1,((rank+1)*(rank+2))/2 920 e2nldd=e2nldd+& 921 & (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+& 922 & gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii)) 923 end do 924 925 ! Terms involving one ucvol derivative (diagonal strain only) 926 ! and one derivative of nonlocal operator 927 ! Loop over strain index 928 do istr2=1,6 929 930 ! rank->rank+2 931 iterm=1 932 temp(:,1:((rank+1)*(rank+2))/2)= & 933 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 934 call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac) 935 e2nl_tmp(istr2)=0.d0 936 do ii=1,((rank+3)*(rank+4))/2 937 e2nl_tmp(istr2)=e2nl_tmp(istr2)-2.d0*& 938 & (dgxds(1,jjs2-1+ii,ia,iproj,isp)*tmpfac(1,ii)+& 939 & dgxds(2,jjs2-1+ii,ia,iproj,isp)*tmpfac(2,ii)) 940 end do 941 ! rank->rank 942 943 iterm=2 944 temp(:,1:((rank+1)*(rank+2))/2)= & 945 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin) 946 call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac) 947 do ii=1,((rank+1)*(rank+2))/2 948 e2nl_tmp(istr2)=e2nl_tmp(istr2)-& 949 & (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+& 950 & gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii)) 951 end do 952 ! DEBUG 953 ! This and subsequent similar debug sections evaluate the 954 ! hermitial conjugate of the contraction immeditely above 955 ! and the comparison was useful for development purposes. 956 ! rank+2->rank 957 ! iterm=3 958 ! temp(:,1:((rank+3)*(rank+4))/2)= & 959 ! dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin) 960 ! call metstr(istr2,rank,iterm,gmet,gprimd,temp,tmpfac) 961 ! e2nl_tmp(istr2)=0.d0 962 ! do ii=1,((rank+1)*(rank+2))/2 963 ! e2nl_tmp(istr2)=e2nl_tmp(istr2)-& 964 ! & (gxa(1,jjs-1+ii,ia,iproj,ispin)*tmpfac(1,ii)+& 965 ! & gxa(2,jjs-1+ii,ia,iproj,ispin)*tmpfac(2,ii)) 966 ! end do 967 ! ENDDEBUG 968 end do 969 970 ! Terms involving two derivatives of the nonlocal operator 971 ! Loop over 2nd strain index 972 do istr2=1,6 973 ! Loop over 1st strain index, upper triangle only 974 do istr1=1,6 975 iest=istr1+(3*natom+6)*(istr2-1) 976 977 ! Accumulate terms corresponding to two derivatives of ucvol 978 ! (simply the nonlocal energy contributin) for both indices 979 ! corresponding to diagonal strains 980 981 if(istr1<=3 .and. istr2<=3) then 982 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nldd 983 end if 984 985 ! Accumulate terms computed above from 1st derivatives 986 ! when one or more indices corresponds to diagonal strain 987 if(istr2<=3) then 988 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl_tmp(istr1) 989 end if 990 if(istr1<=3) then 991 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl_tmp(istr2) 992 end if 993 994 ! rank->rank+4 995 call contstr21(istr2,istr1,rank,gmet,gprimd,e2nl,& 996 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),& 997 & d2gxds2(:,jjs4:jjs4-1+((rank+5)*(rank+6))/2,ia,iproj,isp)) 998 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl 999 1000 ! DEBUG 1001 ! rank+4->rank 1002 ! call contstr22(istr2,istr1,rank,gmet,gprimd,e2nl,& 1003 ! & d2gxds2(:,jjs4:jjs4-1+((rank+5)*(rank+6))/2,ia,iproj,ispin),& 1004 ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)) 1005 ! enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl 1006 ! ENDDEBUG 1007 1008 ! rank->rank+2 1009 call contstr23(istr2,istr1,rank,gmet,gprimd,e2nl,& 1010 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),& 1011 & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp)) 1012 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl 1013 ! DEBUG 1014 1015 ! rank+2->rank 1016 ! call contstr24(istr2,istr1,rank,gmet,gprimd,e2nl,& 1017 ! & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin),& 1018 ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)) 1019 ! enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl 1020 ! ENDDEBUG 1021 1022 ! rank+2->rank+2 1023 if(rank<=2) then 1024 call contstr25(istr2,istr1,rank,gmet,gprimd,e2nl,& 1025 & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp),& 1026 & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp)) 1027 else 1028 call contstr25a(istr2,istr1,rank,gmet,gprimd,e2nl,& 1029 & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp),& 1030 & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp)) 1031 end if 1032 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl 1033 1034 ! rank->rank 1035 call contstr26(istr2,istr1,rank,gmet,gprimd,e2nl,& 1036 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),& 1037 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)) 1038 enlout(iest)= enlout(iest)+wt(ilang,iproj)*e2nl 1039 1040 end do !istr1 1041 1042 ! Contributions to internal strain (one cartesian strain and one 1043 ! reduced-coordinate atomic displacement derivative). 1044 iest=7+3*(ia5-1)+(3*natom+6)*(istr2-1) 1045 1046 ! rank->rank+3 1047 call contistr03(istr2,rank,gmet,gprimd,eisnl,& 1048 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),& 1049 & d2gxdis(:,jjs3:jjs3-1+((rank+4)*(rank+5))/2,ia,iproj,isp)) 1050 enlout(iest:iest+2)= enlout(iest:iest+2)& 1051 & +wt(ilang,iproj)*eisnl(1:3) 1052 1053 ! DEBUG 1054 ! rank+3->rank 1055 ! call contistr30(istr2,rank,gmet,gprimd,eisnl,& 1056 ! & d2gxdis(:,jjs3:jjs3-1+((rank+4)*(rank+5))/2,ia,iproj,ispin),& 1057 ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)) 1058 ! enlout(iest:iest+2)= enlout(iest:iest+2)& 1059 ! & +wt(ilang,iproj)*eisnl(1:3) 1060 ! ENDDEBUG 1061 1062 ! rank+1->rank+2 1063 call contistr12(istr2,rank,gmet,gprimd,eisnl,& 1064 & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,isp),& 1065 & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,isp)) 1066 enlout(iest:iest+2)= enlout(iest:iest+2)& 1067 & +wt(ilang,iproj)*eisnl(1:3) 1068 1069 ! DEBUG 1070 ! rank+2->rank+1 1071 ! call contistr21(istr2,rank,gmet,gprimd,eisnl,& 1072 ! & dgxds(:,jjs2:jjs2-1+((rank+3)*(rank+4))/2,ia,iproj,ispin),& 1073 ! & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,ispin)) 1074 ! enlout(iest:iest+2)= enlout(iest:iest+2)& 1075 ! & +wt(ilang,iproj)*eisnl(1:3) 1076 ! ENDDEBUG 1077 1078 ! rank->rank+1 1079 call contistr01(istr2,rank,gmet,gprimd,eisnl,& 1080 & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin),& 1081 & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,isp)) 1082 enlout(iest:iest+2)= enlout(iest:iest+2)& 1083 & +wt(ilang,iproj)*eisnl(1:3) 1084 ! 1085 ! DEBUG 1086 ! rank+1->rank 1087 ! call contistr10(istr2,rank,gmet,gprimd,eisnl,& 1088 ! & dgxdis(:,jjs1:jjs1-1+((rank+2)*(rank+3))/2,ia,iproj,ispin),& 1089 ! & gxa(:,jjs:jjs-1+((rank+1)*(rank+2))/2,ia,iproj,ispin)) 1090 ! enlout(iest:iest+2)= enlout(iest:iest+2)& 1091 ! & +wt(ilang,iproj)*eisnl(1:3) 1092 ! ENDDEBUG 1093 1094 end do !istr2 1095 end if !choice==6 1096 1097 ! End loop over iproj: 1098 end do 1099 ! End condition of existence of a reference state: 1100 end if 1101 1102 ! End loop over ilang: 1103 end do 1104 1105 ! End loop over ia: 1106 end do 1107 1108 ! Operate with the non-local potential on the projected scalars, 1109 ! in order to get matrix element [NOT needed if only force or stress 1110 ! or dynamical matrix is being computed]. 1111 1112 if (signs==2) then 1113 if(nloalg(2)<=0 .and. choice==2)then 1114 ! Prepare the phase k+q factors for the atoms between ia3 and ia4: 1115 ! they were not prepared previously for nloalg(2)<=0 and choice==2. 1116 call ph1d3d(ia3,ia4,kgout,matblk,natom,npwout,& 1117 & n1,n2,n3,phkxredout,ph1d,ph3dout) 1118 end if 1119 1120 ! call timab(74+choice,1,tsec) 1121 sign=-1 1122 ! The duplication of the opernl calls has the purpose to avoid 1123 ! copying vectout/vectout_s 1124 if(ispin==1.or.nspinor==1)then 1125 if( nloalg(1)==2) then 1126 call opernl2(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,& 1127 & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 1128 & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 1129 & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,& 1130 & ntypat,ph3dout,sign,vectout) 1131 else if( nloalg(1)==3) then 1132 call opernl3(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,& 1133 & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 1134 & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 1135 & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,& 1136 & ntypat,ph3dout,sign,vectout) 1137 else if( nloalg(1)==4) then 1138 call opernl4b(choice,dgxdsfac,dgxdtfac,ffnlout,gmet,gxafac,& 1139 & ia3,idir,indlmn_s,ispinor,itypat,jproj,kgout,kpgout,kptout,& 1140 & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,& 1141 & dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,ntypat,ph3dout,vectout) 1142 end if 1143 else ! if ispin == 2 1144 vectout_s(:,:)=zero 1145 if( nloalg(1)==2) then 1146 call opernl2(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,& 1147 & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 1148 & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 1149 & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,& 1150 & ntypat,ph3dout,sign,vectout_s) 1151 else if( nloalg(1)==3) then 1152 call opernl3(choice,dgxdis,dgxdsfac,d2gxdis,d2gxds2,dgxdtfac,& 1153 & ffnlout,gmet,gxafac,ia3,idir,indlmn_s,ispinor,istwf_k,itypat,& 1154 & jproj,kgout,kpgout,kptout,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 1155 & mlang5,mlang6,mproj,ndgxdt,dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,& 1156 & ntypat,ph3dout,sign,vectout_s) 1157 else if( nloalg(1)==4) then 1158 call opernl4b(choice,dgxds,dgxdtfac,ffnlout,gmet,gxafac,& 1159 & ia3,idir,indlmn_s,ispinor,itypat,jproj,kgout,kpgout,kptout,& 1160 & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,& 1161 & dimffnlout,nincat,nkpgout,nlang,nloalg,npwout,ntypat,ph3dout,vectout_s) 1162 end if 1163 vectout(1,1+npwout:2*npwout)=& 1164 & vectout(1,1+npwout:2*npwout)+vectout_s(1,1:npwout) 1165 vectout(2,1+npwout:2*npwout)=& 1166 & vectout(2,1+npwout:2*npwout)+vectout_s(2,1:npwout) 1167 1168 end if ! end ispin if 1169 end if ! end signs==2 if 1170 1171 ! End loops over spins: 1172 end do 1173 1174 ! End condition of existence of a non-local part for that type of atom: 1175 end if 1176 1177 ! End loop over ispinor: 1178 end do 1179 1180 ! End sum on atom subset loop, over ia3: 1181 end do 1182 1183 ! End atom type loop, over itypat: 1184 ia1=ia2+1 1185 end do 1186 1187 !De-allocate temporary space. 1188 ABI_FREE(ekb_s) 1189 ABI_FREE(gxa) 1190 ABI_FREE(gxafac) 1191 ABI_FREE(dgxds) 1192 ABI_FREE(dgxdt) 1193 ABI_FREE(dgxdtfac) 1194 ABI_FREE(wt) 1195 ABI_FREE(jproj) 1196 ABI_FREE(temp) 1197 ABI_FREE(tmpfac) 1198 ABI_FREE(dgxdsfac) 1199 ABI_FREE(indlmn_s) 1200 !if(choice==6) then 1201 ABI_FREE(dgxdis) 1202 ABI_FREE(d2gxdis) 1203 ABI_FREE(d2gxds2) 1204 !end if 1205 !if(nspinor==2) then 1206 ABI_FREE(dgxds_s) 1207 ABI_FREE(dgxdt_s) 1208 ABI_FREE(gxa_s) 1209 !end if 1210 !if(nspinor==2.and.choice==6) then 1211 ABI_FREE(dgxdis_s) 1212 ABI_FREE(d2gxdis_s) 1213 ABI_FREE(d2gxds2_s) 1214 !end if 1215 if (mpssoang>mpsang) then 1216 ABI_FREE(pauli) 1217 end if 1218 1219 !Restore the original content of the vectin array. 1220 !Note that only the first part was modified 1221 if(istwf_k/=1) then 1222 call scalewf_nonlop(istwf_k,mpi_enreg,npwin,2,vectin) 1223 end if 1224 1225 if (nspinor==2) then 1226 ABI_FREE(vectin_s) 1227 ABI_FREE(vectout_s) 1228 end if 1229 1230 if (mpi_enreg%paral_spinor==1) then 1231 if (size(enlout)>0) call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr) 1232 call xmpi_sum(strsnl,mpi_enreg%comm_spinor,ierr) 1233 call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr) 1234 call xmpi_sum(strsso,mpi_enreg%comm_spinor,ierr) 1235 end if 1236 1237 !Save non-local energy 1238 if((choice==1).and.size(enlout)>0) enlout(1)=enlk ! on test v4/93 size(enlout) can be zero (PMA) 1239 1240 !Do final manipulations to produce strain gradients for 1241 !stress tensor, in cartesian coordinates 1242 if ((choice==3.or.choice==23) .and. signs==1) then 1243 ! Convert strsnl from reduced to cartesian coordinates 1244 strsnl_out(:)=0.d0 1245 call strconv(strsnl,gprimd,strsnl_out) 1246 strsnl(:) = strsnl_out(:) 1247 ! Add diagonal part (fill up first 6 components of enlout with 1248 ! these gradients; elements 7,8,9 of enlout are not used) 1249 enlout(1)=strsnl(1)-enlk 1250 enlout(2)=strsnl(2)-enlk 1251 enlout(3)=strsnl(3)-enlk 1252 enlout(4)=strsnl(4) 1253 enlout(5)=strsnl(5) 1254 enlout(6)=strsnl(6) 1255 ! Eventually, add spin-orbit part due to gradients of 1256 ! antisymmetric tensor 1257 if (mpssoang>mpsang) then 1258 call strsocv(strsso,gprimd,strssoc) 1259 enlout(1:6)=enlout(1:6)+strssoc(:) 1260 end if 1261 end if 1262 1263 !DEBUG 1264 !write(std_out,*)' nonlop_pl: exit ' 1265 !ENDDEBUG 1266 1267 contains
ABINIT/scalewf_nonlop [ Functions ]
NAME
scalewf_nonlop
FUNCTION
At the start of nonlop (or similar routines), as well as its end, the wavefunctions, when stored with istwfk/=2, need to be scaled (by a factor of 2 or 1/2), except for the G=0 component. Only the first spinor component is to be modified.
INPUTS
istwf_k=storage mode of the vector mpi_enreg=information about MPI parallelization npw=number of planewaves option=1 multiply by 2 =2 multiply by 1/2
OUTPUT
(see side effects)
SIDE EFFECTS
vect(2,npw)=vector that is rescaled
NOTES
XG030513 : MPIWF One should pay attention to the G=0 component, that will be only one one proc...
SOURCE
1431 subroutine scalewf_nonlop(istwf_k,mpi_enreg,npw,option,vect) 1432 1433 !Arguments ------------------------------------ 1434 !This type is defined in defs_mpi 1435 !scalars 1436 integer,intent(in) :: istwf_k,npw,option 1437 type(MPI_type),intent(in) :: mpi_enreg 1438 !arrays 1439 real(dp),intent(inout) :: vect(2,npw) 1440 1441 !Local variables------------------------------- 1442 !scalars 1443 integer :: ipw 1444 real(dp) :: scale 1445 character(len=500) :: msg 1446 1447 ! ************************************************************************* 1448 1449 DBG_ENTER("COLL") 1450 1451 if(istwf_k/=1)then 1452 1453 if(option/=1 .and. option/=2)then 1454 write(msg,'(a,a,a,i0)')& 1455 'The argument option should be 1 or 2,',ch10,& 1456 'however, option=',option 1457 ABI_BUG(msg) 1458 end if 1459 1460 scale=two 1461 if(option==2)scale=half 1462 1463 ! Storage for the Gamma point. The component of the G=0 vector 1464 ! should not be scaled, and no G=0 imaginary part is allowed. 1465 if(istwf_k==2)then 1466 if (mpi_enreg%me_g0==1) then 1467 vect(2,1)=zero 1468 !$OMP PARALLEL DO 1469 do ipw=2,npw 1470 vect(1,ipw)=scale*vect(1,ipw) 1471 vect(2,ipw)=scale*vect(2,ipw) 1472 end do 1473 !$OMP END PARALLEL DO 1474 else 1475 !$OMP PARALLEL DO 1476 do ipw=1,npw 1477 vect(1,ipw)=scale*vect(1,ipw) 1478 vect(2,ipw)=scale*vect(2,ipw) 1479 end do 1480 !$OMP END PARALLEL DO 1481 end if 1482 end if 1483 1484 ! Other storage modes, for k points with time-reversal symmetry. 1485 ! All components should be scaled. 1486 if(istwf_k>2)then 1487 !$OMP PARALLEL DO 1488 do ipw=1,npw 1489 vect(1,ipw)=scale*vect(1,ipw) 1490 vect(2,ipw)=scale*vect(2,ipw) 1491 end do 1492 !$OMP END PARALLEL DO 1493 end if 1494 1495 end if ! istwf_k/=1 1496 1497 DBG_EXIT("COLL") 1498 1499 end subroutine scalewf_nonlop
ABINIT/strsocv [ Functions ]
NAME
strsocv
FUNCTION
Convert from antisymmetric storage mode 3x3x3 rank3 tensor in reduced coordinates "red" to symmetric storage mode 3x3 rank2 tensor in cartesian coordinates "cart", using metric tensor "gprimd".
INPUTS
red(6,3)=3x3x3 tensor in antisymmetric storage mode, reduced coordinates gprimd(3,3)=reciprocal space dimensional primitive translations
OUTPUT
cart(6)=3x3 tensor in symmetric storage mode, cartesian coordinates
NOTES
This routine is used to compute spin-orbit stress tensor. red is antisymmetric in first two indices and stored as 11 22 33 32 31 21. cart is symmetric and stored as 11 22 33 32 31 21. cart(1,1) & = & & red(i,j,2) G(3,i) G(1,j) + red(i,j,3) G(1,i) G(2,j) \nonumber cart(2,2) & = & & red(i,j,1) G(2,i) G(3,j) + red(i,j,3) G(1,i) G(2,j) \nonumber cart(3,3) & = & & red(i,j,1) G(2,i) G(3,j) + red(i,j,2) G(3,i) G(1,j) \nonumber cart(3,2) & = & 0.5 ( & red(i,j,3) G(1,i) G(3,j) + red(i,j,2) G(2,i) G(1,j)) \nonumber cart(3,1) & = & 0.5 ( & red(i,j,3) G(3,i) G(2,j) + red(i,j,1) G(2,i) G(1,j)) \nonumber cart(2,1) & = & 0.5 ( & red(i,j,2) G(3,i) G(2,j) + red(i,j,1) G(1,i) G(3,j)) \end{eqnarray} }}
SOURCE
1353 subroutine strsocv(red,gprimd,cart) 1354 1355 !Arguments ------------------------------------ 1356 !arrays 1357 real(dp),intent(in) :: gprimd(3,3),red(6,3) 1358 real(dp),intent(out) :: cart(6) 1359 1360 !Local variables------------------------------- 1361 !scalars 1362 integer :: ii,jj 1363 !arrays 1364 real(dp) :: work(3,3,3) 1365 1366 ! ************************************************************************* 1367 1368 do ii=1,3 1369 work(1,1,ii)=0.d0 1370 work(2,2,ii)=0.d0 1371 work(3,3,ii)=0.d0 1372 work(3,2,ii)=red(4,ii) ; work(2,3,ii)=-red(4,ii) 1373 work(3,1,ii)=red(5,ii) ; work(1,3,ii)=-red(5,ii) 1374 work(2,1,ii)=red(6,ii) ; work(1,2,ii)=-red(6,ii) 1375 end do 1376 1377 cart(:)=0.d0 1378 do jj=1,3 1379 do ii=1,3 1380 cart(1)=cart(1)+work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj)& 1381 & +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj) 1382 cart(2)=cart(2)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)& 1383 & +work(ii,jj,3)*gprimd(1,ii)*gprimd(2,jj) 1384 cart(3)=cart(3)+work(ii,jj,1)*gprimd(2,ii)*gprimd(3,jj)& 1385 & +work(ii,jj,2)*gprimd(3,ii)*gprimd(1,jj) 1386 cart(4)=cart(4)+work(ii,jj,3)*gprimd(1,ii)*gprimd(3,jj)& 1387 & +work(ii,jj,2)*gprimd(2,ii)*gprimd(1,jj) 1388 cart(5)=cart(5)+work(ii,jj,3)*gprimd(3,ii)*gprimd(2,jj)& 1389 & +work(ii,jj,1)*gprimd(2,ii)*gprimd(1,jj) 1390 cart(6)=cart(6)+work(ii,jj,2)*gprimd(3,ii)*gprimd(2,jj)& 1391 & +work(ii,jj,1)*gprimd(1,ii)*gprimd(3,jj) 1392 end do 1393 end do 1394 cart(4)=0.5d0*cart(4) 1395 cart(5)=0.5d0*cart(5) 1396 cart(6)=0.5d0*cart(6) 1397 1398 end subroutine strsocv
ABINIT/trace2 [ Functions ]
NAME
trace2
FUNCTION
Sum indices to compute trace of rank 2 tensor gxa related to l=2 $trace=sum_{i,j} {gxa(i,j) gmet(i,j)}$
INPUTS
gxa(2,6) = $sum_{G} {e^(2 \pi i G cdot t} {{f_2(|k+G|)} \over {|k+G|^2}} (k+G) cdot (k+G) C(G_{nk})}$ gmet(3,3)=(symmetric) metric tensor in reciprocal space (bohr^-2)
OUTPUT
trace(2)=sum_{i,j} {gxa(i,j) gmet(i,j)}$ (Re and Im).
NOTES
Here index 6 refers to vector components of (k+G) but note tensor is symmetric=>only 6 components. The components are given in the order 11 22 33 32 31 21. The initial 2 handles the Re and Im parts.
SOURCE
1293 subroutine trace2(gxa,gmet,trace) 1294 1295 !Arguments ------------------------------------ 1296 !arrays 1297 real(dp),intent(in) :: gmet(3,3),gxa(2,6) 1298 real(dp),intent(out) :: trace(2) 1299 1300 !Local variables------------------------------- 1301 !scalars 1302 integer :: reim 1303 1304 ! ************************************************************************* 1305 1306 !Write out index summation, Re and Im parts 1307 do reim=1,2 1308 trace(reim)=gxa(reim,1)*gmet(1,1)+gxa(reim,2)*gmet(2,2)+& 1309 & gxa(reim,3)*gmet(3,3)+& 1310 & 2.0d0*(gxa(reim,4)*gmet(3,2)+gxa(reim,5)*gmet(3,1)+& 1311 & gxa(reim,6)*gmet(2,1)) 1312 end do 1313 1314 end subroutine trace2