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