TABLE OF CONTENTS
ABINIT/m_opernlb_ylm [ Modules ]
NAME
m_opernlb_ylm
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (MT) 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_opernlb_ylm 22 23 use defs_basis 24 use m_abicore 25 use m_errors 26 #if defined HAVE_OPENMP 27 use OMP_LIB 28 #endif 29 30 implicit none 31 32 private
ABINIT/opernlb_ylm [ Functions ]
NAME
opernlb_ylm
FUNCTION
* Operate with the non-local part of the hamiltonian, from projected scalars to reciprocal space. * Operate with the non-local projectors and the overlap matrix, from projected scalars to reciprocal space.
INPUTS
choice=chooses possible output (see below)
cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
2 if <p_lmn|c> scalars are complex
cplex_dgxdt(ndgxdt_fac) = used only when cplex = 1
cplex_dgxdt(i)=1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Vnl (NL operator)
dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfacrelated to Sij (overlap)
dimffnl=second dimension of ffnl
ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))= reduced projected scalars related to Sij (overlap)
ia3=gives the number of the first atom in the subset presently treated
idir=direction of the - atom to be moved in the case (choice=2,signs=2) or (choice=22,signs=2)
- k point direction in the case (choice=5, 51, 52 and signs=2)
- strain component (1:6) in the case (choice=2,signs=2) or (choice=6,signs=1)
- strain component (1:9) in the case (choice=33,signs=2)
- (1:9) components to specify the atom to be moved and the second q-gradient
direction in the case (choice=25,signs=2)
indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
kpg(npw,nkpg)=(k+G) components (if nkpg=3).
(k+G) Cartesian components for choice=33
matblk=dimension of the array ph3d
ndgxdtfac=second dimension of dgxdtfac
nincat=number of atoms in the subset here treated
nkpg=second dimension of array kpg (0 or 3)
nlmn=number of (l,m,n) numbers for current type of atom
nloalg(3)=governs the choice of the algorithm for non-local operator.
npw=number of plane waves in reciprocal space
nspinor=number of spinorial components of the wavefunctions (on current proc)
paw_opt= define the nonlocal operator concerned with:
paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
paw_opt=3 : PAW overlap matrix (Sij)
paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
ph3d(2,npw,matblk)=three-dimensional phase factors
[qdir]= optional, direction of the q-gradient (only for choice=22, choice=25 and choice=33)
ucvol=unit cell volume (bohr^3)
OUTPUT
(see side effects)
SIDE EFFECTS
--if (paw_opt=0)
vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
or one of its derivatives to the input vect.
if (choice=22) <G|d2V_nonlocal/d(atm. pos)dq|vect_in> (at q=0)
if (choice=25) <G|d3V_nonlocal/d(atm. pos)dqdq|vect_in> (at q=0)
if (choice=33) <G|d2V_nonlocal/d(strain)dq|vect_in> (at q=0)
--if (paw_opt=0, 1 or 4)
vect(2,npwout*nspinor)=result of the aplication of the concerned operator
or one of its derivatives to the input vect.:
if (choice=1) <G|V_nonlocal|vect_in>
if (choice=2) <G|dV_nonlocal/d(atm. pos)|vect_in>
if (choice=3) <G|dV_nonlocal/d(strain)|vect_in>
if (choice=5) <G|dV_nonlocal/d(k)|vect_in>
if (choice=51) <G|d(right)V_nonlocal/d(k)|vect_in>
if (choice=52) <G|d(left)V_nonlocal/d(k)|vect_in>
if (choice=53) <G|d(twist)V_nonlocal/d(k)|vect_in>
if (choice=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|vect_in>
if (choice=8) <G|d2V_nonlocal/d(k)d(k)|vect_in>
if (choice=81) <G|d[d(right)V_nonlocal/d(k)]/d(k)|vect_in>
if (paw_opt=2)
vect(2,npwout*nspinor)=final vector in reciprocal space:
if (choice=1) <G|V_nonlocal-lamdba.(I+S)|vect_in> (note: not including <G|I|c>)
if (choice=2) <G|d[V_nonlocal-lamdba.(I+S)]/d(atm. pos)|vect_in>
if (choice=3) <G|d[V_nonlocal-lamdba.(I+S)]/d(strain)|vect_in>
if (choice=5) <G|d[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=51) <G|d(right)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=52) <G|d(left)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=53) <G|d(twist)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|vect_in>
if (choice=8) <G|d2[V_nonlocal-lamdba.(I+S)]/d(k)d(k)|vect_in>
if (choice=81) <G|d[d(right[V_nonlocal-lamdba.(I+S)]/d(k)]/d(k)|vect_in>
--if (paw_opt=3 or 4)
svect(2,npwout*nspinor)=result of the aplication of Sij (overlap matrix)
or one of its derivatives to the input vect.:
if (choice=1) <G|I+S|vect_in> (note: not including <G|I|c>)
if (choice=2) <G|dS/d(atm. pos)|vect_in>
if (choice=3) <G|dS/d(strain)|vect_in>
if (choice=5) <G|dS/d(k)|vect_in>
if (choice=51) <G|d(right)S/d(k)|vect_in>
if (choice=52) <G|d(left)S/d(k)|vect_in>
if (choice=53) <G|d(twist)S/d(k)|vect_in>
if (choice=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|vect_in>
if (choice=7) <G|sum_i[p_i><p_i]|vect_in>
if (choice=8) <G|d2S/d(k)d(k)|vect_in>
if (choice=81) <G|d[d(right)S/d(k)]/d(k)|vect_in>
NOTES
1-The openMP version is different from the standard version: the standard version is more effifient on one CPU core. 2-Operate for one type of atom, and within this given type of atom, for a subset of at most nincat atoms.
SOURCE
152 subroutine opernlb_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,& 153 & d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnl,ffnl,gxfac,gxfac_sij,& 154 & ia3,idir,indlmn,kpg,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpg,nlmn,nloalg,npw,& 155 & nspinor,paw_opt,ph3d,svect,ucvol,vect,qdir) 156 157 !Arguments ------------------------------------ 158 !scalars 159 integer,intent(in) :: choice,cplex,cplex_fac,dimffnl,ia3,idir,matblk,ndgxdtfac,nd2gxdtfac,nincat 160 integer,intent(in) :: nkpg,nlmn,npw,nspinor,paw_opt 161 integer,intent(in),optional :: qdir 162 real(dp),intent(in) :: ucvol 163 !arrays 164 integer,intent(in) :: cplex_dgxdt(ndgxdtfac),cplex_d2gxdt(nd2gxdtfac),indlmn(6,nlmn),nloalg(3) 165 real(dp),intent(in) :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor) 166 real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor) 167 real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat*(paw_opt/3),nspinor) 168 real(dp),intent(in) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat*(paw_opt/3),nspinor) 169 real(dp),intent(in) :: ffnl(npw,dimffnl,nlmn),gxfac(cplex_fac,nlmn,nincat,nspinor) 170 real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3)) 171 real(dp),intent(in) :: kpg(npw,nkpg),ph3d(2,npw,matblk) 172 real(dp),intent(inout) :: svect(:,:),vect(:,:) 173 !Local variables------------------------------- 174 !Arrays 175 !scalars 176 integer :: fdb,fdf,ia,ialpha,iaph3d,ibeta,ic,idelta,idelgam,igamma 177 integer :: ii,il,ilmn,ipw,ipwshft,ispinor,jc,nthreads,ffnl_dir1,ffnl_dir(3) 178 real(dp) :: scale,two_piinv,wt 179 logical :: parity 180 !arrays 181 integer,parameter :: ffnl_dir_dat(6)=(/3,4,4,2,2,3/) 182 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 183 integer,parameter :: idir1(9)=(/1,1,1,2,2,2,3,3,3/),idir2(9)=(/1,2,3,1,2,3,1,2,3/) 184 integer,parameter :: nalpha(9)=(/1,2,3,3,3,2,2,1,1/),nbeta(9)=(/1,2,3,2,1,1,3,3,2/) 185 real(dp),allocatable :: d2gxdtfac_(:,:,:),d2gxdtfacs_(:,:,:),dgxdtfac_(:,:,:),dgxdtfacs_(:,:,:),gxfac_(:,:),gxfacs_(:,:) 186 ! real(dp),allocatable :: kpg(:,:) 187 complex(dpc),allocatable :: ztab(:) 188 189 ! ************************************************************************* 190 191 DBG_ENTER("COLL") 192 193 !Nothing to do when choice=4, 6 or 23 194 if (choice==4.or.choice==6.or.choice==23) return 195 196 !DDK not compatible with istwkf > 1 197 if(cplex==1.and.(any(cplex_dgxdt(:)==2).or.any(cplex_d2gxdt(:)==2)))then 198 ABI_BUG("opernlb_ylm+ddk not compatible with istwfk>1") 199 end if 200 201 !Inits 202 wt=four_pi/sqrt(ucvol) 203 nthreads=1 204 #if defined HAVE_OPENMP 205 nthreads=OMP_GET_NUM_THREADS() 206 #endif 207 208 if (paw_opt/=3) then 209 ABI_MALLOC(gxfac_,(2,nlmn)) 210 gxfac_(:,:)=zero 211 if (choice>1) then 212 ABI_MALLOC(dgxdtfac_,(2,ndgxdtfac,nlmn)) 213 if(ndgxdtfac>0) dgxdtfac_(:,:,:)=zero 214 end if 215 if (choice==54.or.choice==8.or.choice==81.or.choice==33) then 216 ABI_MALLOC(d2gxdtfac_,(2,nd2gxdtfac,nlmn)) 217 if(nd2gxdtfac>0) d2gxdtfac_(:,:,:)=zero 218 end if 219 end if 220 if (paw_opt>=3) then 221 ABI_MALLOC(gxfacs_,(2,nlmn)) 222 gxfacs_(:,:)=zero 223 if (choice>1) then 224 ABI_MALLOC(dgxdtfacs_,(2,ndgxdtfac,nlmn)) 225 if (ndgxdtfac>0) dgxdtfacs_(:,:,:)=zero 226 end if 227 if (choice==54.or.choice==8.or.choice==81) then 228 ABI_MALLOC(d2gxdtfacs_,(2,nd2gxdtfac,nlmn)) 229 if (nd2gxdtfac>0) d2gxdtfacs_(:,:,:)=zero 230 end if 231 end if 232 233 if (choice==33) two_piinv=1.0_dp/two_pi 234 235 if (opernlb_counter>=0) then 236 opernlb_counter = opernlb_counter + 1 237 if (paw_opt==4) opernlb_counter = opernlb_counter + 1 238 end if 239 240 ABI_MALLOC(ztab,(npw)) 241 242 !========================================================================== 243 !========== STANDARD VERSION ============================================== 244 !========================================================================== 245 if (nthreads==1) then 246 247 ! Loop on spinorial components 248 do ispinor=1,nspinor 249 ipwshft=(ispinor-1)*npw 250 251 ! Loop on atoms (blocking) 252 do ia=1,nincat 253 iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1 254 ! Scale gxfac with 4pi/sqr(omega).(-i)^l 255 if (paw_opt/=3) then 256 do ilmn=1,nlmn 257 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 258 scale=wt;if (il>1) scale=-scale 259 if (parity) then 260 gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor) 261 if (cplex_fac==1) gxfac_(2,ilmn)=zero 262 else 263 gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor) 264 if (cplex_fac==2) then 265 gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor) 266 else 267 gxfac_(1,ilmn)=zero 268 end if 269 end if! parity 270 end do ! ilmn 271 if (choice>1) then 272 do ilmn=1,nlmn 273 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 274 scale=wt;if (il>1) scale=-scale 275 if (parity) then 276 if(cplex_fac==2)then 277 dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor) 278 else 279 do ii=1,ndgxdtfac 280 ic = cplex_dgxdt(ii) ; jc = 3-ic 281 dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 282 dgxdtfac_(jc,ii,ilmn)=zero 283 end do 284 end if 285 else 286 if(cplex_fac==2)then 287 do ii=1,ndgxdtfac 288 dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor) 289 dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 290 end do 291 else 292 do ii=1,ndgxdtfac 293 ic = cplex_dgxdt(ii) ; jc = 3-ic 294 dgxdtfac_(ic,ii,ilmn)=zero 295 if(ic==1)then 296 dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 297 else 298 dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 299 end if 300 end do 301 end if 302 end if 303 end do 304 end if ! choice>1 305 if (choice==54.or.choice==8.or.choice==81.or.choice==33) then 306 do ilmn=1,nlmn 307 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 308 scale=wt;if (il>1) scale=-scale 309 if (parity) then 310 if(cplex_fac==2)then 311 d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor) 312 else 313 do ii=1,nd2gxdtfac 314 ic = cplex_d2gxdt(ii) ; jc = 3-ic 315 d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 316 d2gxdtfac_(jc,ii,ilmn)=zero 317 end do 318 end if 319 else 320 if(cplex_fac==2)then 321 do ii=1,nd2gxdtfac 322 d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor) 323 d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 324 end do 325 else 326 do ii=1,nd2gxdtfac 327 ic = cplex_d2gxdt(ii) ; jc = 3-ic 328 d2gxdtfac_(ic,ii,ilmn)=zero 329 if(ic==1)then 330 d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 331 else 332 d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 333 end if 334 end do 335 end if 336 end if 337 end do! ilmn 338 end if ! choice 54 or 8 or 81 or 33 339 end if ! paw_opt /= 3 340 341 ! Scale gxfac_sij with 4pi/sqr(omega).(-i)^l 342 if (paw_opt>=3) then 343 344 do ilmn=1,nlmn 345 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 346 scale=wt;if (il>1) scale=-scale 347 if (parity) then 348 gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor) 349 if (cplex==1) gxfacs_(2,ilmn)=zero 350 else 351 gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor) 352 if (cplex==2) then 353 gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor) 354 else 355 gxfacs_(1,ilmn)=zero 356 end if 357 end if 358 end do 359 if (choice>1) then 360 do ilmn=1,nlmn 361 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 362 scale=wt;if (il>1) scale=-scale 363 if (parity) then 364 if(cplex==2)then 365 dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn) = scale * dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor) 366 else 367 do ii=1,ndgxdtfac 368 ic = cplex_dgxdt(ii) ; jc = 3-ic 369 dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 370 dgxdtfacs_(jc,ii,ilmn)=zero 371 end do 372 end if 373 else 374 if(cplex==2)then 375 do ii=1,ndgxdtfac 376 dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor) 377 dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 378 end do 379 else 380 do ii=1,ndgxdtfac 381 ic = cplex_dgxdt(ii) ; jc = 3-ic 382 dgxdtfacs_(ic,ii,ilmn)=zero 383 if(ic==1)then 384 dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 385 else 386 dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 387 end if 388 end do 389 end if 390 end if 391 end do ! ilmn 392 end if ! choice>1 393 if (choice==54.or.choice==8.or.choice==81) then 394 do ilmn=1,nlmn 395 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 396 scale=wt;if (il>1) scale=-scale 397 if (parity) then 398 if(cplex==2)then 399 d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor) 400 else 401 do ii=1,nd2gxdtfac 402 ic = cplex_d2gxdt(ii) ; jc = 3-ic 403 d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 404 d2gxdtfacs_(jc,ii,ilmn)=zero 405 end do 406 end if 407 else 408 if(cplex==2)then 409 do ii=1,nd2gxdtfac 410 d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor) 411 d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 412 end do 413 else 414 do ii=1,nd2gxdtfac 415 ic = cplex_d2gxdt(ii) ; jc = 3-ic 416 d2gxdtfacs_(ic,ii,ilmn)=zero 417 if(ic==1)then 418 d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 419 else 420 d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 421 end if 422 end do 423 end if 424 end if ! parity 425 end do ! ilmn 426 end if ! choice == 54 or 8 or 81 427 end if ! paw_opt >= 3 428 429 ! Compute <g|Vnl|c> (or derivatives) for each plane wave: 430 431 if (paw_opt/=3) then 432 433 ztab(:)=czero 434 435 ! ------ 436 if (choice==1) then ! <g|Vnl|c> 437 do ilmn=1,nlmn 438 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 439 end do 440 end if 441 442 ! ------ 443 if (choice==2) then ! derivative w.r.t. atm. pos 444 do ilmn=1,nlmn 445 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp) 446 end do 447 ztab(:)=two_pi*kpg(:,idir)*ztab(:) 448 do ilmn=1,nlmn 449 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 450 end do 451 end if 452 453 ! ------ 454 if (choice==22) then ! mixed derivative w.r.t. atm. pos and q vector (at q=0) 455 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+qdir 456 if (idir==qdir) then 457 do ilmn=1,nlmn 458 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 459 end do 460 end if 461 do ilmn=1,nlmn 462 ztab(:)=ztab(:)+kpg(:,idir)*ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 463 ztab(:)=ztab(:)-ffnl(:,ffnl_dir1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 464 end do 465 ztab(:)=ztab(:)*two_pi 466 end if 467 468 ! ------ 469 if (choice==25) then ! mixed derivative w.r.t. atm. pos and two q vectors (at q=0) 470 !Use same notation as the notes for clarity 471 ialpha=nalpha(idir) 472 idelta=nbeta(idir) 473 igamma=qdir 474 idelgam=gamma(idelta,igamma) 475 if (ialpha==igamma) then 476 do ilmn=1,nlmn 477 ztab(:)=ztab(:)+ffnl(:,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 478 end do 479 end if 480 if (ialpha==idelta) then 481 do ilmn=1,nlmn 482 ztab(:)=ztab(:)+ffnl(:,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 483 end do 484 end if 485 do ilmn=1,nlmn 486 ztab(:)=ztab(:)+kpg(:,ialpha)*ffnl(:,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 487 ztab(:)=ztab(:)-ffnl(:,4+idelgam,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 488 end do 489 ztab(:)=ztab(:)*two_pi 490 end if 491 492 ! ------ 493 if (choice==3) then ! derivative w.r.t. strain 494 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 495 if (idir<=3) then 496 do ilmn=1,nlmn 497 ztab(:)=ztab(:)+ffnl(:,1,ilmn)& 498 & *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp)& 499 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 500 end do 501 else 502 do ilmn=1,nlmn 503 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)& 504 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 505 end do 506 end if 507 end if 508 509 ! ------ 510 if (choice==33) then ! mixed derivative w.r.t. strain and q vector (at q=0) 511 !Use same notation as the notes for clarity 512 ibeta=nalpha(idir) 513 idelta=nbeta(idir) 514 igamma=qdir 515 idelgam=gamma(idelta,igamma) 516 if (ibeta==igamma) then 517 do ilmn=1,nlmn 518 ztab(:)=ztab(:)+onehalf*ffnl(:,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 519 ztab(:)=ztab(:)+half*ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) 520 end do 521 end if 522 if (ibeta==idelta) then 523 do ilmn=1,nlmn 524 ztab(:)=ztab(:)+onehalf*ffnl(:,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 525 ztab(:)=ztab(:)+half*ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 526 end do 527 end if 528 do ilmn=1,nlmn 529 ztab(:)=ztab(:)+kpg(:,ibeta)*ffnl(:,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 530 ztab(:)=ztab(:)+ffnl(:,1+idelta,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp) 531 ztab(:)=ztab(:)+ffnl(:,1+igamma,ilmn)*cmplx(d2gxdtfac_(1,2,ilmn),d2gxdtfac_(2,2,ilmn),kind=dp) 532 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfac_(1,3,ilmn),d2gxdtfac_(2,3,ilmn),kind=dp) 533 end do 534 ztab(:)=ztab(:)*two_piinv 535 end if 536 537 ! ------ 538 if (choice==5) then ! full derivative w.r.t. k 539 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 540 do ilmn=1,nlmn 541 ztab(:)=ztab(:)+ffnl(:,1 ,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)& 542 & +ffnl(:,ffnl_dir1,ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2 ,ilmn),kind=dp) 543 end do 544 end if 545 546 ! ------ 547 if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi> 548 do ilmn=1,nlmn 549 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 550 end do 551 end if 552 553 ! ------ 554 if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi> 555 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 556 do ilmn=1,nlmn 557 ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 558 end do 559 end if 560 561 ! ------ 562 if (choice==53) then ! twist derivative: <G|dp_i/dk_(idir+1)>V_ij<dp_j/dk_(idir+2)|psi> 563 fdf = ffnl_dir_dat(2*idir-1) 564 fdb = ffnl_dir_dat(2*idir) 565 do ilmn=1,nlmn 566 ztab(:)=ztab(:) + & 567 & ffnl(:,fdf,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) 568 end do 569 end if 570 571 ! ------ 572 if (choice==54) then ! mixed derivative w.r.t. atm. pos and (right) k 573 do ilmn=1,nlmn 574 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp) 575 end do 576 ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:) 577 do ilmn=1,nlmn 578 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp) 579 end do 580 end if 581 582 ! ------ 583 if (choice==8) then ! full second order derivative w.r.t. k 584 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 585 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 586 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 587 do ilmn=1,nlmn 588 ztab(:)=ztab(:) & 589 & +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2, ilmn),kind=dp)& 590 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 591 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)& 592 & +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 593 end do 594 end if 595 596 ! ------ 597 if (choice==81) then 598 ! partial second order derivative w.r.t. k 599 ! full derivative w.r.t. k1, right derivative w.r.t. k2 600 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 601 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 602 do ilmn=1,nlmn 603 ztab(:)=ztab(:) & 604 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 605 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 606 end do 607 end if 608 609 ! ------ 610 ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp) 611 612 vect(1,1+ipwshft:npw+ipwshft)=vect(1,1+ipwshft:npw+ipwshft)+real(ztab(:)) 613 vect(2,1+ipwshft:npw+ipwshft)=vect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:)) 614 615 end if ! paw_opt /= 3 616 617 ! Compute <g|S|c> (or derivatives) for each plane wave: 618 619 if (paw_opt>=3) then 620 621 ztab(:)=czero 622 623 ! ------ 624 if (choice==1) then ! <g|S|c> 625 do ilmn=1,nlmn 626 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 627 end do 628 end if 629 630 ! ------ 631 if (choice==2) then ! derivative w.r.t. atm. pos 632 do ilmn=1,nlmn 633 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp) 634 end do 635 ztab(:)=two_pi*kpg(:,idir)*ztab(:) 636 do ilmn=1,nlmn 637 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 638 end do 639 end if 640 641 ! ------ 642 if (choice==3) then ! derivative w.r.t. strain 643 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 644 if (idir<=3) then 645 do ilmn=1,nlmn 646 ztab(:)=ztab(:)+ffnl(:,1,ilmn)& 647 & *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)& 648 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 649 end do 650 else 651 do ilmn=1,nlmn 652 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)& 653 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 654 end do 655 end if 656 end if 657 658 ! ------ 659 if (choice==5) then ! full derivative w.r.t. k 660 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 661 do ilmn=1,nlmn 662 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)& 663 & +ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 664 end do 665 end if 666 667 ! ------ 668 if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi> 669 do ilmn=1,nlmn 670 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 671 end do 672 end if 673 674 ! ------ 675 if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi> 676 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 677 do ilmn=1,nlmn 678 ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 679 end do 680 end if 681 682 ! ------ 683 if (choice==53) then ! twist derivative: <G|dp_i/dk_(idir+1)>S_ij<dp_j/dk_(idir+2)|psi> 684 fdf = ffnl_dir_dat(2*idir-1) 685 fdb = ffnl_dir_dat(2*idir) 686 do ilmn=1,nlmn 687 ztab(:)=ztab(:) + & 688 & ffnl(:,fdf,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp) 689 end do 690 end if 691 692 ! ------ 693 if (choice==54) then ! mixed derivative w.r.t. atm. pos and k 694 do ilmn=1,nlmn 695 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp) 696 end do 697 ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:) 698 do ilmn=1,nlmn 699 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp) 700 end do 701 end if 702 703 ! ------ 704 if (choice==8) then ! full second order derivative w.r.t. k 705 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 706 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 707 do ilmn=1,nlmn 708 ztab(:)=ztab(:) & 709 & +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx( gxfacs_(1, ilmn), gxfacs_(2, ilmn),kind=dp)& 710 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 711 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)& 712 & +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 713 end do 714 end if 715 716 ! ------ 717 if (choice==81) then 718 ! partial second order derivative w.r.t. k 719 ! full derivative w.r.t. k1, right derivative w.r.t. k2 720 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 721 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 722 do ilmn=1,nlmn 723 ztab(:)=ztab(:) & 724 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 725 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 726 end do 727 end if 728 729 730 ! ------ 731 ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp) 732 svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:)) 733 svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:)) 734 end if ! paw_opt >= 3 735 736 ! End loop on atoms 737 end do 738 end do ! End loop on spinors 739 740 741 ! ========================================================================== 742 ! ========== OPENMP VERSION ================================================ 743 ! ========================================================================== 744 else 745 746 ! Loop on spinorial components 747 do ispinor=1,nspinor 748 ipwshft=(ispinor-1)*npw 749 750 ! Loop on atoms (blocking) 751 do ia=1,nincat 752 iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1 753 754 ! Scale gxfac with 4pi/sqr(omega).(-i)^l 755 if (paw_opt/=3) then 756 !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc) 757 !$OMP DO 758 do ilmn=1,nlmn 759 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 760 scale=wt;if (il>1) scale=-scale 761 if (parity) then 762 gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor) 763 if (cplex_fac==1) gxfac_(2,ilmn)=zero 764 else 765 gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor) 766 if (cplex_fac==2) then 767 gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor) 768 else 769 gxfac_(1,ilmn)=zero 770 end if 771 end if 772 end do 773 !$OMP END DO 774 if (choice>1) then 775 !$OMP DO 776 do ilmn=1,nlmn 777 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 778 scale=wt;if (il>1) scale=-scale 779 if (parity) then 780 if(cplex_fac==2)then 781 dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor) 782 else 783 do ii=1,ndgxdtfac 784 ic = cplex_dgxdt(ii) ; jc = 3-ic 785 dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 786 dgxdtfac_(jc,ii,ilmn)=zero 787 end do 788 end if 789 else 790 if(cplex_fac==2)then 791 do ii=1,ndgxdtfac 792 dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 793 dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor) 794 end do 795 else 796 do ii=1,ndgxdtfac 797 ic = cplex_dgxdt(ii) ; jc = 3-ic 798 dgxdtfac_(ic,ii,ilmn)=zero 799 if(ic==1)then 800 dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 801 else 802 dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 803 end if 804 end do 805 end if 806 end if 807 end do 808 !$OMP END DO 809 end if 810 if (choice==54.or.choice==8.or.choice==81.or.choice==33) then 811 !$OMP DO 812 do ilmn=1,nlmn 813 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 814 scale=wt;if (il>1) scale=-scale 815 if (parity) then 816 if(cplex_fac==2)then 817 d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor) 818 else 819 do ii=1,nd2gxdtfac 820 ic = cplex_d2gxdt(ii) ; jc = 3-ic 821 d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 822 d2gxdtfac_(jc,ii,ilmn)=zero 823 end do 824 end if 825 else 826 if(cplex_fac==2)then 827 do ii=1,nd2gxdtfac 828 d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor) 829 d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 830 end do 831 else 832 do ii=1,nd2gxdtfac 833 ic = cplex_d2gxdt(ii) ; jc = 3-ic 834 d2gxdtfac_(ic,ii,ilmn)=zero 835 if(ic==1)then 836 d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 837 else 838 d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 839 end if 840 end do 841 end if 842 end if 843 end do 844 !$OMP END DO 845 end if 846 !$OMP END PARALLEL 847 end if 848 849 ! Scale gxfac_sij with 4pi/sqr(omega).(-i)^l 850 if (paw_opt>=3) then 851 !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc) 852 !$OMP DO 853 do ilmn=1,nlmn 854 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 855 scale=wt;if (il>1) scale=-scale 856 if (parity) then 857 gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor) 858 if (cplex==1) gxfacs_(2,ilmn)=zero 859 else 860 gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor) 861 if (cplex==2) then 862 gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor) 863 else 864 gxfacs_(1,ilmn)=zero 865 end if 866 end if 867 end do 868 !$OMP END DO 869 if (choice>1) then 870 !$OMP DO 871 do ilmn=1,nlmn 872 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 873 scale=wt;if (il>1) scale=-scale 874 if (parity) then 875 if(cplex==2)then 876 dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn)=scale*dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor) 877 else 878 do ii=1,ndgxdtfac 879 ic = cplex_dgxdt(ii) ; jc = 3-ic 880 dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 881 dgxdtfacs_(jc,ii,ilmn)=zero 882 end do 883 end if 884 else 885 if(cplex==2)then 886 do ii=1,ndgxdtfac 887 dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 888 dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor) 889 end do 890 else 891 do ii=1,ndgxdtfac 892 ic = cplex_dgxdt(ii) ; jc = 3-ic 893 dgxdtfacs_(ic,ii,ilmn)=zero 894 if(ic==1)then 895 dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 896 else 897 dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 898 end if 899 end do 900 end if 901 end if 902 end do 903 !$OMP END DO 904 end if 905 if (choice==54.or.choice==8.or.choice==81) then 906 !$OMP DO 907 do ilmn=1,nlmn 908 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 909 scale=wt;if (il>1) scale=-scale 910 if (parity) then 911 if(cplex==2)then 912 d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor) 913 else 914 do ii=1,nd2gxdtfac 915 ic = cplex_d2gxdt(ii) ; jc = 3-ic 916 d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 917 d2gxdtfacs_(jc,ii,ilmn)=zero 918 end do 919 end if 920 else 921 if(cplex==2)then 922 do ii=1,nd2gxdtfac 923 d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor) 924 d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 925 end do 926 else 927 do ii=1,nd2gxdtfac 928 ic = cplex_d2gxdt(ii) ; jc = 3-ic 929 d2gxdtfacs_(ic,ii,ilmn)=zero 930 if(ic==1)then 931 d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 932 else 933 d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 934 end if 935 end do 936 end if 937 end if 938 end do 939 !$OMP END DO 940 end if 941 !$OMP END PARALLEL 942 end if 943 944 ! Compute <g|Vnl|c> (or derivatives) for each plane wave: 945 if (paw_opt/=3) then 946 !$OMP PARALLEL PRIVATE(ipw,ilmn,fdf,fdb,ffnl_dir1) 947 948 ! ------ 949 if (choice==1) then ! <g|Vnl|c> 950 !$OMP DO 951 do ipw=1,npw 952 ztab(ipw)=czero 953 do ilmn=1,nlmn 954 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 955 end do 956 end do 957 !$OMP END DO 958 959 ! ------ 960 else if (choice==2) then ! derivative w.r.t. atm. pos 961 !$OMP DO 962 do ipw=1,npw 963 ztab(ipw)=czero 964 do ilmn=1,nlmn 965 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp) 966 end do 967 ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw) 968 do ilmn=1,nlmn 969 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 970 end do 971 end do 972 !$OMP END DO 973 974 ! ------ 975 else if (choice==22) then ! mixed derivative w.r.t. atm. pos and q vector (at q=0) 976 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+qdir 977 !$OMP DO 978 do ipw=1,npw 979 ztab(ipw)=czero 980 if (idir==qdir) then 981 do ilmn=1,nlmn 982 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 983 end do 984 end if 985 do ilmn=1,nlmn 986 ztab(ipw)=ztab(ipw)+kpg(ipw,idir)*ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 987 ztab(ipw)=ztab(ipw)-ffnl(ipw,ffnl_dir1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 988 end do 989 ztab(ipw)=ztab(ipw)*two_pi 990 end do 991 !$OMP END DO 992 993 ! ------ 994 else if (choice==25) then ! mixed derivative w.r.t. atm. pos and thwo q vectors (at q=0) 995 !Use same notation as the notes for clarity 996 ialpha=nalpha(idir) 997 idelta=nbeta(idir) 998 igamma=qdir 999 idelgam=gamma(idelta,igamma) 1000 !$OMP DO 1001 do ipw=1,npw 1002 ztab(ipw)=czero 1003 if (ialpha==igamma) then 1004 do ilmn=1,nlmn 1005 ztab(ipw)=ztab(ipw)+ffnl(ipw,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1006 end do 1007 end if 1008 if (ialpha==idelta) then 1009 do ilmn=1,nlmn 1010 ztab(ipw)=ztab(ipw)+ffnl(ipw,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1011 end do 1012 end if 1013 do ilmn=1,nlmn 1014 ztab(ipw)=ztab(ipw)+kpg(ipw,ialpha)*ffnl(ipw,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1015 ztab(ipw)=ztab(ipw)-ffnl(ipw,4+idelgam,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 1016 end do 1017 ztab(ipw)=ztab(ipw)*two_pi 1018 end do 1019 !$OMP END DO 1020 ! ------ 1021 else if (choice==3) then ! derivative w.r.t. strain 1022 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1023 if (idir<=3) then 1024 !$OMP DO 1025 do ipw=1,npw 1026 ztab(ipw)=czero 1027 do ilmn=1,nlmn 1028 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) & 1029 & *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp) & 1030 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1031 end do 1032 end do 1033 !$OMP END DO 1034 else 1035 !$OMP DO 1036 do ipw=1,npw 1037 ztab(ipw)=czero 1038 do ilmn=1,nlmn 1039 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) & 1040 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1041 end do 1042 end do 1043 !$OMP END DO 1044 end if 1045 1046 ! ------ 1047 else if (choice==33) then ! mixed derivative w.r.t. strain and q vector (at q=0) 1048 !Use same notation as the notes for clarity 1049 ibeta=nalpha(idir) 1050 idelta=nbeta(idir) 1051 igamma=qdir 1052 idelgam=gamma(idelta,igamma) 1053 !$OMP DO 1054 do ipw=1,npw 1055 ztab(ipw)=czero 1056 if (ibeta==igamma) then 1057 do ilmn=1,nlmn 1058 ztab(ipw)=ztab(ipw)+onehalf*ffnl(ipw,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1059 ztab(ipw)=ztab(ipw)+half*ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) 1060 end do 1061 end if 1062 if (ibeta==idelta) then 1063 do ilmn=1,nlmn 1064 ztab(ipw)=ztab(ipw)+onehalf*ffnl(ipw,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1065 ztab(ipw)=ztab(ipw)+half*ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 1066 end do 1067 end if 1068 do ilmn=1,nlmn 1069 ztab(ipw)=ztab(ipw)+kpg(ipw,ibeta)*ffnl(ipw,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1070 ztab(ipw)=ztab(ipw)+ffnl(ipw,1+idelta,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp) 1071 ztab(ipw)=ztab(ipw)+ffnl(ipw,1+igamma,ilmn)*cmplx(d2gxdtfac_(1,2,ilmn),d2gxdtfac_(2,2,ilmn),kind=dp) 1072 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfac_(1,3,ilmn),d2gxdtfac_(2,3,ilmn),kind=dp) 1073 end do 1074 ztab(ipw)=ztab(ipw)*two_piinv 1075 end do 1076 !$OMP END DO 1077 1078 ! ------ 1079 else if (choice==5) then ! full derivative w.r.t. k 1080 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1081 !$OMP DO 1082 do ipw=1,npw 1083 ztab(ipw)=czero 1084 do ilmn=1,nlmn 1085 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) & 1086 & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1087 end do 1088 end do 1089 !$OMP END DO 1090 1091 ! ------ 1092 else if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi> 1093 !$OMP DO 1094 do ipw=1,npw 1095 ztab(ipw)=czero 1096 do ilmn=1,nlmn 1097 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 1098 end do 1099 end do 1100 !$OMP END DO 1101 1102 ! ------ 1103 else if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi> 1104 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1105 !$OMP DO 1106 do ipw=1,npw 1107 ztab(ipw)=czero 1108 do ilmn=1,nlmn 1109 ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 1110 end do 1111 end do 1112 !$OMP END DO 1113 1114 ! ------ 1115 else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>V<dp/dk_(idir+2)|psi> 1116 fdf = ffnl_dir_dat(2*idir-1) 1117 fdb = ffnl_dir_dat(2*idir) 1118 !$OMP DO 1119 do ipw=1,npw 1120 ztab(ipw)=czero 1121 do ilmn=1,nlmn 1122 ztab(ipw)=ztab(ipw) & 1123 & +ffnl(ipw,fdf,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) 1124 end do 1125 end do 1126 !$OMP END DO 1127 1128 ! ------ 1129 else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k 1130 !$OMP DO 1131 do ipw=1,npw 1132 ztab(ipw)=czero 1133 do ilmn=1,nlmn 1134 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp) 1135 end do 1136 ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw) 1137 do ilmn=1,nlmn 1138 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp) 1139 end do 1140 end do 1141 !$OMP END DO 1142 1143 ! ------ 1144 else if (choice==8) then ! full second order derivative w.r.t. k 1145 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1146 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 1147 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 1148 !$OMP DO 1149 do ipw=1,npw 1150 ztab(ipw)=czero 1151 do ilmn=1,nlmn 1152 ztab(ipw)=ztab(ipw) & 1153 & +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2, ilmn),kind=dp)& 1154 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 1155 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)& 1156 & +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 1157 end do 1158 end do 1159 !$OMP END DO 1160 1161 ! ------ 1162 else if (choice==81) then 1163 ! partial second order derivative w.r.t. k 1164 ! full derivative w.r.t. k1, right derivative w.r.t. k2 1165 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1166 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 1167 !$OMP DO 1168 do ipw=1,npw 1169 ztab(ipw)=czero 1170 do ilmn=1,nlmn 1171 ztab(ipw)=ztab(ipw) & 1172 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 1173 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 1174 end do 1175 end do 1176 !$OMP END DO 1177 1178 ! ------ 1179 else 1180 !$OMP WORKSHARE 1181 ztab(:)=czero 1182 !$OMP END WORKSHARE 1183 end if 1184 1185 1186 ! ------ 1187 !$OMP DO 1188 do ipw=1,npw 1189 ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp) 1190 vect(1,ipw+ipwshft)=vect(1,ipw+ipwshft)+real(ztab(ipw)) 1191 vect(2,ipw+ipwshft)=vect(2,ipw+ipwshft)+aimag(ztab(ipw)) 1192 end do 1193 !$OMP END DO 1194 1195 !$OMP END PARALLEL 1196 end if 1197 1198 ! Compute <g|S|c> (or derivatives) for each plane wave: 1199 if (paw_opt>=3) then 1200 !$OMP PARALLEL PRIVATE(ilmn,ipw,fdf,fdb,ffnl_dir1) 1201 1202 ! ------ 1203 if (choice==1) then ! <g|S|c> 1204 !$OMP DO 1205 do ipw=1,npw 1206 ztab(ipw)=czero 1207 do ilmn=1,nlmn 1208 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1209 end do 1210 end do 1211 !$OMP END DO 1212 1213 ! ------ 1214 else if (choice==2) then ! derivative w.r.t. atm. pos 1215 !$OMP DO 1216 do ipw=1,npw 1217 ztab(ipw)=czero 1218 do ilmn=1,nlmn 1219 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp) 1220 end do 1221 ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw) 1222 do ilmn=1,nlmn 1223 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 1224 end do 1225 end do 1226 !$OMP END DO 1227 1228 ! ------ 1229 else if (choice==3) then ! derivative w.r.t. strain 1230 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1231 if (idir<=3) then 1232 !$OMP DO 1233 do ipw=1,npw 1234 ztab(ipw)=czero 1235 do ilmn=1,nlmn 1236 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) & 1237 & *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)& 1238 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1239 end do 1240 end do 1241 !$OMP END DO 1242 else 1243 !$OMP DO 1244 do ipw=1,npw 1245 ztab(ipw)=czero 1246 do ilmn=1,nlmn 1247 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) & 1248 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1249 end do 1250 end do 1251 !$OMP END DO 1252 end if 1253 1254 ! ------ 1255 else if (choice==5) then ! full derivative w.r.t. k 1256 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1257 !$OMP DO 1258 do ipw=1,npw 1259 ztab(ipw)=czero 1260 do ilmn=1,nlmn 1261 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) & 1262 & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1263 end do 1264 end do 1265 !$OMP END DO 1266 1267 ! ------ 1268 else if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi> 1269 !$OMP DO 1270 do ipw=1,npw 1271 ztab(ipw)=czero 1272 do ilmn=1,nlmn 1273 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 1274 end do 1275 end do 1276 !$OMP END DO 1277 1278 ! ------ 1279 else if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi> 1280 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1281 !$OMP DO 1282 do ipw=1,npw 1283 ztab(ipw)=czero 1284 do ilmn=1,nlmn 1285 ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1286 end do 1287 end do 1288 !$OMP END DO 1289 1290 ! ------ 1291 else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>S<dp/dk_(idir+2)|psi> 1292 fdf = ffnl_dir_dat(2*idir-1) 1293 fdb = ffnl_dir_dat(2*idir) 1294 !$OMP DO 1295 do ipw=1,npw 1296 ztab(ipw)=czero 1297 do ilmn=1,nlmn 1298 ztab(ipw)=ztab(ipw) & 1299 & +ffnl(ipw,fdf,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp) 1300 end do 1301 end do 1302 !$OMP END DO 1303 1304 ! ------ 1305 else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k 1306 !$OMP DO 1307 do ipw=1,npw 1308 ztab(ipw)=czero 1309 do ilmn=1,nlmn 1310 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp) 1311 end do 1312 ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw) 1313 do ilmn=1,nlmn 1314 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp) 1315 end do 1316 end do 1317 !$OMP END DO 1318 1319 ! ------ 1320 else if (choice==8) then ! full second order derivative w.r.t. k 1321 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1322 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 1323 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 1324 !$OMP DO 1325 do ipw=1,npw 1326 ztab(ipw)=czero 1327 do ilmn=1,nlmn 1328 ztab(ipw)=ztab(ipw) & 1329 & +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx( gxfacs_(1, ilmn), gxfacs_(2, ilmn),kind=dp)& 1330 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 1331 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)& 1332 & +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 1333 end do 1334 end do 1335 !$OMP END DO 1336 1337 ! ------ 1338 else if (choice==81) then 1339 ! partial second order derivative w.r.t. k 1340 ! full derivative w.r.t. k1, right derivative w.r.t. k2 1341 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1342 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 1343 !$OMP DO 1344 do ipw=1,npw 1345 ztab(ipw)=czero 1346 do ilmn=1,nlmn 1347 ztab(ipw)=ztab(ipw) & 1348 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 1349 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 1350 end do 1351 end do 1352 !$OMP END DO 1353 1354 ! ------ 1355 else 1356 !$OMP WORKSHARE 1357 ztab(:)=czero 1358 !$OMP END WORKSHARE 1359 end if 1360 1361 1362 ! ------ 1363 ! The OMP WORKSHARE directive doesn't have a good performance with Intel Compiler 1364 ! !$OMP WORKSHARE 1365 ! ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp) 1366 ! !$OMP END WORKSHARE 1367 ! !$OMP WORKSHARE 1368 ! svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:)) 1369 ! svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:)) 1370 ! !$OMP END WORKSHARE 1371 !$OMP DO 1372 do ipw=1,npw 1373 ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp) 1374 svect(1,ipw+ipwshft)=svect(1,ipw+ipwshft)+real(ztab(ipw)) 1375 svect(2,ipw+ipwshft)=svect(2,ipw+ipwshft)+aimag(ztab(ipw)) 1376 end do 1377 !$OMP END DO 1378 !$OMP END PARALLEL 1379 end if 1380 1381 ! End loop on atoms 1382 end do 1383 ! End loop on spinors 1384 end do 1385 1386 ! ========================================================================== 1387 end if 1388 1389 ABI_FREE(ztab) 1390 1391 if (paw_opt/=3) then 1392 ABI_FREE(gxfac_) 1393 if (choice>1) then 1394 ABI_FREE(dgxdtfac_) 1395 end if 1396 if (choice==54.or.choice==8.or.choice==81.or.choice==33) then 1397 ABI_FREE(d2gxdtfac_) 1398 end if 1399 end if 1400 if (paw_opt>=3) then 1401 ABI_FREE(gxfacs_) 1402 if (choice>1) then 1403 ABI_FREE(dgxdtfacs_) 1404 end if 1405 if (choice==54.or.choice==8.or.choice==81) then 1406 ABI_FREE(d2gxdtfacs_) 1407 end if 1408 end if 1409 1410 DBG_EXIT("COLL") 1411 1412 #if !defined HAVE_OPENMP 1413 !Fake use of unused variable 1414 if (.false.) write(std_out,*) ipw 1415 #endif 1416 1417 end subroutine opernlb_ylm