TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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),
- 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)
indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
kpg(npw,nkpg)=(k+G) components (if nkpg=3)
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
ucvol=unit cell volume (bohr^3)
OUTPUT
(see side effects)
SIDE EFFECTS
--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.
PARENTS
nonlop_ylm
CHILDREN
SOURCE
113 #if defined HAVE_CONFIG_H 114 #include "config.h" 115 #endif 116 117 #include "abi_common.h" 118 119 subroutine opernlb_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,& 120 & d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnl,ffnl,gxfac,gxfac_sij,& 121 & ia3,idir,indlmn,kpg,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpg,nlmn,nloalg,npw,& 122 & nspinor,paw_opt,ph3d,svect,ucvol,vect) 123 124 use defs_basis 125 use m_profiling_abi 126 use m_errors 127 #if defined HAVE_OPENMP 128 use OMP_LIB 129 #endif 130 131 !This section has been created automatically by the script Abilint (TD). 132 !Do not modify the following lines by hand. 133 #undef ABI_FUNC 134 #define ABI_FUNC 'opernlb_ylm' 135 !End of the abilint section 136 137 implicit none 138 139 !Arguments ------------------------------------ 140 !scalars 141 142 integer,intent(in) :: choice,cplex,cplex_fac,dimffnl,ia3,idir,matblk,ndgxdtfac,nd2gxdtfac,nincat 143 integer,intent(in) :: nkpg,nlmn,npw,nspinor,paw_opt 144 real(dp),intent(in) :: ucvol 145 !arrays 146 integer,intent(in) :: cplex_dgxdt(ndgxdtfac),cplex_d2gxdt(nd2gxdtfac),indlmn(6,nlmn),nloalg(3) 147 real(dp),intent(in) :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor) 148 real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor) 149 real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat*(paw_opt/3),nspinor) 150 real(dp),intent(in) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat*(paw_opt/3),nspinor) 151 real(dp),intent(in) :: ffnl(npw,dimffnl,nlmn),gxfac(cplex_fac,nlmn,nincat,nspinor) 152 real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3)) 153 real(dp),intent(in) :: kpg(npw,nkpg),ph3d(2,npw,matblk) 154 real(dp),intent(inout) :: svect(:,:),vect(:,:) 155 !Local variables------------------------------- 156 !Arrays 157 !scalars 158 integer :: fdb,fdf,ia,iaph3d,ic,ii,il,ilmn,ipw,ipwshft,ispinor,jc,nthreads,ffnl_dir1,ffnl_dir(3) 159 real(dp) :: scale,wt 160 logical :: parity 161 !arrays 162 integer,parameter :: ffnl_dir_dat(6)=(/3,4,4,2,2,3/) 163 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 164 integer,parameter :: idir1(9)=(/1,1,1,2,2,2,3,3,3/),idir2(9)=(/1,2,3,1,2,3,1,2,3/) 165 real(dp),allocatable :: d2gxdtfac_(:,:,:),d2gxdtfacs_(:,:,:),dgxdtfac_(:,:,:),dgxdtfacs_(:,:,:),gxfac_(:,:),gxfacs_(:,:) 166 complex(dpc),allocatable :: ztab(:) 167 168 ! ************************************************************************* 169 170 DBG_ENTER("COLL") 171 172 !Nothing to do when choice=4, 6 or 23 173 if (choice==4.or.choice==6.or.choice==23) return 174 175 !DDK not compatible with istwkf > 1 176 if(cplex==1.and.(any(cplex_dgxdt(:)==2).or.any(cplex_d2gxdt(:)==2)))then 177 MSG_BUG("opernlb_ylm+ddk not compatible with istwfk>1") 178 end if 179 180 !Inits 181 wt=four_pi/sqrt(ucvol) 182 nthreads=1 183 #if defined HAVE_OPENMP 184 nthreads=OMP_GET_NUM_THREADS() 185 #endif 186 187 if (paw_opt/=3) then 188 ABI_ALLOCATE(gxfac_,(2,nlmn)) 189 gxfac_(:,:)=zero 190 if (choice>1) then 191 ABI_ALLOCATE(dgxdtfac_,(2,ndgxdtfac,nlmn)) 192 if(ndgxdtfac>0) dgxdtfac_(:,:,:)=zero 193 end if 194 if (choice==54.or.choice==8.or.choice==81) then 195 ABI_ALLOCATE(d2gxdtfac_,(2,nd2gxdtfac,nlmn)) 196 if(nd2gxdtfac>0) d2gxdtfac_(:,:,:)=zero 197 end if 198 end if 199 if (paw_opt>=3) then 200 ABI_ALLOCATE(gxfacs_,(2,nlmn)) 201 gxfacs_(:,:)=zero 202 if (choice>1) then 203 ABI_ALLOCATE(dgxdtfacs_,(2,ndgxdtfac,nlmn)) 204 if (ndgxdtfac>0) dgxdtfacs_(:,:,:)=zero 205 end if 206 if (choice==54.or.choice==8.or.choice==81) then 207 ABI_ALLOCATE(d2gxdtfacs_,(2,nd2gxdtfac,nlmn)) 208 if (nd2gxdtfac>0) d2gxdtfacs_(:,:,:)=zero 209 end if 210 end if 211 212 ABI_ALLOCATE(ztab,(npw)) 213 214 !========================================================================== 215 !========== STANDARD VERSION ============================================== 216 !========================================================================== 217 if (nthreads==1) then 218 219 ! Loop on spinorial components 220 do ispinor=1,nspinor 221 ipwshft=(ispinor-1)*npw 222 223 ! Loop on atoms (blocking) 224 do ia=1,nincat 225 iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1 226 ! Scale gxfac with 4pi/sqr(omega).(-i)^l 227 if (paw_opt/=3) then 228 do ilmn=1,nlmn 229 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 230 scale=wt;if (il>1) scale=-scale 231 if (parity) then 232 gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor) 233 if (cplex_fac==1) gxfac_(2,ilmn)=zero 234 else 235 gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor) 236 if (cplex_fac==2) then 237 gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor) 238 else 239 gxfac_(1,ilmn)=zero 240 end if 241 end if 242 end do 243 if (choice>1) then 244 do ilmn=1,nlmn 245 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 246 scale=wt;if (il>1) scale=-scale 247 if (parity) then 248 if(cplex_fac==2)then 249 dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor) 250 else 251 do ii=1,ndgxdtfac 252 ic = cplex_dgxdt(ii) ; jc = 3-ic 253 dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 254 dgxdtfac_(jc,ii,ilmn)=zero 255 end do 256 end if 257 else 258 if(cplex_fac==2)then 259 do ii=1,ndgxdtfac 260 dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor) 261 dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 262 end do 263 else 264 do ii=1,ndgxdtfac 265 ic = cplex_dgxdt(ii) ; jc = 3-ic 266 dgxdtfac_(ic,ii,ilmn)=zero 267 if(ic==1)then 268 dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 269 else 270 dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 271 end if 272 end do 273 end if 274 end if 275 end do 276 end if 277 if (choice==54.or.choice==8.or.choice==81) then 278 do ilmn=1,nlmn 279 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 280 scale=wt;if (il>1) scale=-scale 281 if (parity) then 282 if(cplex_fac==2)then 283 d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor) 284 else 285 do ii=1,nd2gxdtfac 286 ic = cplex_d2gxdt(ii) ; jc = 3-ic 287 d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 288 d2gxdtfac_(jc,ii,ilmn)=zero 289 end do 290 end if 291 else 292 if(cplex_fac==2)then 293 do ii=1,nd2gxdtfac 294 d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor) 295 d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 296 end do 297 else 298 do ii=1,nd2gxdtfac 299 ic = cplex_d2gxdt(ii) ; jc = 3-ic 300 d2gxdtfac_(ic,ii,ilmn)=zero 301 if(ic==1)then 302 d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 303 else 304 d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 305 end if 306 end do 307 end if 308 end if 309 end do 310 end if 311 end if 312 313 ! Scale gxfac_sij with 4pi/sqr(omega).(-i)^l 314 if (paw_opt>=3) then 315 do ilmn=1,nlmn 316 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 317 scale=wt;if (il>1) scale=-scale 318 if (parity) then 319 gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor) 320 if (cplex==1) gxfacs_(2,ilmn)=zero 321 else 322 gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor) 323 if (cplex==2) then 324 gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor) 325 else 326 gxfacs_(1,ilmn)=zero 327 end if 328 end if 329 end do 330 if (choice>1) then 331 do ilmn=1,nlmn 332 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 333 scale=wt;if (il>1) scale=-scale 334 if (parity) then 335 if(cplex==2)then 336 dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn)=scale*dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor) 337 else 338 do ii=1,ndgxdtfac 339 ic = cplex_dgxdt(ii) ; jc = 3-ic 340 dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 341 dgxdtfacs_(jc,ii,ilmn)=zero 342 end do 343 end if 344 else 345 if(cplex==2)then 346 do ii=1,ndgxdtfac 347 dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor) 348 dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 349 end do 350 else 351 do ii=1,ndgxdtfac 352 ic = cplex_dgxdt(ii) ; jc = 3-ic 353 dgxdtfacs_(ic,ii,ilmn)=zero 354 if(ic==1)then 355 dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 356 else 357 dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 358 end if 359 end do 360 end if 361 end if 362 end do 363 end if 364 if (choice==54.or.choice==8.or.choice==81) then 365 do ilmn=1,nlmn 366 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 367 scale=wt;if (il>1) scale=-scale 368 if (parity) then 369 if(cplex==2)then 370 d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor) 371 else 372 do ii=1,nd2gxdtfac 373 ic = cplex_d2gxdt(ii) ; jc = 3-ic 374 d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 375 d2gxdtfacs_(jc,ii,ilmn)=zero 376 end do 377 end if 378 else 379 if(cplex==2)then 380 do ii=1,nd2gxdtfac 381 d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor) 382 d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 383 end do 384 else 385 do ii=1,nd2gxdtfac 386 ic = cplex_d2gxdt(ii) ; jc = 3-ic 387 d2gxdtfacs_(ic,ii,ilmn)=zero 388 if(ic==1)then 389 d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 390 else 391 d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 392 end if 393 end do 394 end if 395 end if 396 end do 397 end if 398 end if 399 400 ! Compute <g|Vnl|c> (or derivatives) for each plane wave: 401 402 if (paw_opt/=3) then 403 404 ztab(:)=czero 405 406 ! ------ 407 if (choice==1) then ! <g|Vnl|c> 408 do ilmn=1,nlmn 409 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 410 end do 411 end if 412 413 ! ------ 414 if (choice==2) then ! derivative w.r.t. atm. pos 415 do ilmn=1,nlmn 416 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp) 417 end do 418 ztab(:)=two_pi*kpg(:,idir)*ztab(:) 419 do ilmn=1,nlmn 420 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 421 end do 422 end if 423 424 ! ------ 425 if (choice==3) then ! derivative w.r.t. strain 426 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 427 if (idir<=3) then 428 do ilmn=1,nlmn 429 ztab(:)=ztab(:)+ffnl(:,1,ilmn)& 430 & *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp)& 431 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 432 end do 433 else 434 do ilmn=1,nlmn 435 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)& 436 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 437 end do 438 end if 439 end if 440 441 ! ------ 442 if (choice==5) then ! full derivative w.r.t. k 443 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 444 do ilmn=1,nlmn 445 ztab(:)=ztab(:)+ffnl(:,1 ,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)& 446 & +ffnl(:,ffnl_dir1,ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2 ,ilmn),kind=dp) 447 end do 448 end if 449 450 ! ------ 451 if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi> 452 do ilmn=1,nlmn 453 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 454 end do 455 end if 456 457 ! ------ 458 if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi> 459 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 460 do ilmn=1,nlmn 461 ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 462 end do 463 end if 464 465 ! ------ 466 if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>V<dp/dk_(idir-1)|psi> 467 ! -<G|dp/dk_(idir-1)>V<dp/dk_(idir+1)|psi> 468 fdf = ffnl_dir_dat(2*idir-1) 469 fdb = ffnl_dir_dat(2*idir) 470 do ilmn=1,nlmn 471 ztab(:)=ztab(:) + & 472 & ffnl(:,fdf,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) - & 473 & ffnl(:,fdb,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 474 end do 475 end if 476 477 ! ------ 478 if (choice==54) then ! mixed derivative w.r.t. atm. pos and (right) k 479 do ilmn=1,nlmn 480 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp) 481 end do 482 ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:) 483 do ilmn=1,nlmn 484 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp) 485 end do 486 end if 487 488 ! ------ 489 if (choice==8) then ! full second order derivative w.r.t. k 490 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 491 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 492 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 493 do ilmn=1,nlmn 494 ztab(:)=ztab(:) & 495 & +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2, ilmn),kind=dp)& 496 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 497 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)& 498 & +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 499 end do 500 end if 501 502 ! ------ 503 if (choice==81) then 504 ! partial second order derivative w.r.t. k 505 ! full derivative w.r.t. k1, right derivative w.r.t. k2 506 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 507 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 508 do ilmn=1,nlmn 509 ztab(:)=ztab(:) & 510 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 511 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 512 end do 513 end if 514 515 ! ------ 516 ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp) 517 518 vect(1,1+ipwshft:npw+ipwshft)=vect(1,1+ipwshft:npw+ipwshft)+real(ztab(:)) 519 vect(2,1+ipwshft:npw+ipwshft)=vect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:)) 520 521 end if 522 523 ! Compute <g|S|c> (or derivatives) for each plane wave: 524 525 if (paw_opt>=3) then 526 527 ztab(:)=czero 528 529 ! ------ 530 if (choice==1) then ! <g|S|c> 531 do ilmn=1,nlmn 532 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 533 end do 534 end if 535 536 ! ------ 537 if (choice==2) then ! derivative w.r.t. atm. pos 538 do ilmn=1,nlmn 539 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp) 540 end do 541 ztab(:)=two_pi*kpg(:,idir)*ztab(:) 542 do ilmn=1,nlmn 543 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 544 end do 545 end if 546 547 ! ------ 548 if (choice==3) then ! derivative w.r.t. strain 549 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 550 if (idir<=3) then 551 do ilmn=1,nlmn 552 ztab(:)=ztab(:)+ffnl(:,1,ilmn)& 553 & *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)& 554 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 555 end do 556 else 557 do ilmn=1,nlmn 558 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)& 559 & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 560 end do 561 end if 562 end if 563 564 ! ------ 565 if (choice==5) then ! full derivative w.r.t. k 566 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 567 do ilmn=1,nlmn 568 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)& 569 & +ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 570 end do 571 end if 572 573 ! ------ 574 if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi> 575 do ilmn=1,nlmn 576 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 577 end do 578 end if 579 580 ! ------ 581 if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi> 582 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 583 do ilmn=1,nlmn 584 ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 585 end do 586 end if 587 588 ! ------ 589 if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>S<dp/dk_(idir-1)|psi> 590 ! -<G|dp/dk_(idir-1)>S<dp/dk_(idir+1)|psi> 591 fdf = ffnl_dir_dat(2*idir-1) 592 fdb = ffnl_dir_dat(2*idir) 593 do ilmn=1,nlmn 594 ztab(:)=ztab(:) + & 595 & ffnl(:,fdf,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp) - & 596 & ffnl(:,fdb,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 597 end do 598 end if 599 600 ! ------ 601 if (choice==54) then ! mixed derivative w.r.t. atm. pos and k 602 do ilmn=1,nlmn 603 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp) 604 end do 605 ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:) 606 do ilmn=1,nlmn 607 ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp) 608 end do 609 end if 610 611 ! ------ 612 if (choice==8) then ! full second order derivative w.r.t. k 613 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 614 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 615 do ilmn=1,nlmn 616 ztab(:)=ztab(:) & 617 & +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx( gxfacs_(1, ilmn), gxfacs_(2, ilmn),kind=dp)& 618 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 619 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)& 620 & +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 621 end do 622 end if 623 624 ! ------ 625 if (choice==81) then 626 ! partial second order derivative w.r.t. k 627 ! full derivative w.r.t. k1, right derivative w.r.t. k2 628 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 629 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 630 do ilmn=1,nlmn 631 ztab(:)=ztab(:) & 632 & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 633 & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 634 end do 635 end if 636 637 638 ! ------ 639 ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp) 640 svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:)) 641 svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:)) 642 end if 643 644 ! End loop on atoms 645 end do 646 end do ! End loop on spinors 647 648 649 ! ========================================================================== 650 ! ========== OPENMP VERSION ================================================ 651 ! ========================================================================== 652 else 653 654 ! Loop on spinorial components 655 do ispinor=1,nspinor 656 ipwshft=(ispinor-1)*npw 657 658 ! Loop on atoms (blocking) 659 do ia=1,nincat 660 iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1 661 662 ! Scale gxfac with 4pi/sqr(omega).(-i)^l 663 if (paw_opt/=3) then 664 !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc) 665 !$OMP DO 666 do ilmn=1,nlmn 667 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 668 scale=wt;if (il>1) scale=-scale 669 if (parity) then 670 gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor) 671 if (cplex_fac==1) gxfac_(2,ilmn)=zero 672 else 673 gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor) 674 if (cplex_fac==2) then 675 gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor) 676 else 677 gxfac_(1,ilmn)=zero 678 end if 679 end if 680 end do 681 !$OMP END DO 682 if (choice>1) then 683 !$OMP DO 684 do ilmn=1,nlmn 685 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 686 scale=wt;if (il>1) scale=-scale 687 if (parity) then 688 if(cplex_fac==2)then 689 dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor) 690 else 691 do ii=1,ndgxdtfac 692 ic = cplex_dgxdt(ii) ; jc = 3-ic 693 dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 694 dgxdtfac_(jc,ii,ilmn)=zero 695 end do 696 end if 697 else 698 if(cplex_fac==2)then 699 do ii=1,ndgxdtfac 700 dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 701 dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor) 702 end do 703 else 704 do ii=1,ndgxdtfac 705 ic = cplex_dgxdt(ii) ; jc = 3-ic 706 dgxdtfac_(ic,ii,ilmn)=zero 707 if(ic==1)then 708 dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 709 else 710 dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor) 711 end if 712 end do 713 end if 714 end if 715 end do 716 !$OMP END DO 717 end if 718 if (choice==54.or.choice==8.or.choice==81) then 719 !$OMP DO 720 do ilmn=1,nlmn 721 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 722 scale=wt;if (il>1) scale=-scale 723 if (parity) then 724 if(cplex_fac==2)then 725 d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor) 726 else 727 do ii=1,nd2gxdtfac 728 ic = cplex_d2gxdt(ii) ; jc = 3-ic 729 d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 730 d2gxdtfac_(jc,ii,ilmn)=zero 731 end do 732 end if 733 else 734 if(cplex_fac==2)then 735 do ii=1,nd2gxdtfac 736 d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor) 737 d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 738 end do 739 else 740 do ii=1,nd2gxdtfac 741 ic = cplex_d2gxdt(ii) ; jc = 3-ic 742 d2gxdtfac_(ic,ii,ilmn)=zero 743 if(ic==1)then 744 d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 745 else 746 d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor) 747 end if 748 end do 749 end if 750 end if 751 end do 752 !$OMP END DO 753 end if 754 !$OMP END PARALLEL 755 end if 756 757 ! Scale gxfac_sij with 4pi/sqr(omega).(-i)^l 758 if (paw_opt>=3) then 759 !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc) 760 !$OMP DO 761 do ilmn=1,nlmn 762 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 763 scale=wt;if (il>1) scale=-scale 764 if (parity) then 765 gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor) 766 if (cplex==1) gxfacs_(2,ilmn)=zero 767 else 768 gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor) 769 if (cplex==2) then 770 gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor) 771 else 772 gxfacs_(1,ilmn)=zero 773 end if 774 end if 775 end do 776 !$OMP END DO 777 if (choice>1) then 778 !$OMP DO 779 do ilmn=1,nlmn 780 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 781 scale=wt;if (il>1) scale=-scale 782 if (parity) then 783 if(cplex==2)then 784 dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn)=scale*dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor) 785 else 786 do ii=1,ndgxdtfac 787 ic = cplex_dgxdt(ii) ; jc = 3-ic 788 dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 789 dgxdtfacs_(jc,ii,ilmn)=zero 790 end do 791 end if 792 else 793 if(cplex==2)then 794 do ii=1,ndgxdtfac 795 dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 796 dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor) 797 end do 798 else 799 do ii=1,ndgxdtfac 800 ic = cplex_dgxdt(ii) ; jc = 3-ic 801 dgxdtfacs_(ic,ii,ilmn)=zero 802 if(ic==1)then 803 dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 804 else 805 dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor) 806 end if 807 end do 808 end if 809 end if 810 end do 811 !$OMP END DO 812 end if 813 if (choice==54.or.choice==8.or.choice==81) then 814 !$OMP DO 815 do ilmn=1,nlmn 816 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 817 scale=wt;if (il>1) scale=-scale 818 if (parity) then 819 if(cplex==2)then 820 d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor) 821 else 822 do ii=1,nd2gxdtfac 823 ic = cplex_d2gxdt(ii) ; jc = 3-ic 824 d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 825 d2gxdtfacs_(jc,ii,ilmn)=zero 826 end do 827 end if 828 else 829 if(cplex==2)then 830 do ii=1,nd2gxdtfac 831 d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor) 832 d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 833 end do 834 else 835 do ii=1,nd2gxdtfac 836 ic = cplex_d2gxdt(ii) ; jc = 3-ic 837 d2gxdtfacs_(ic,ii,ilmn)=zero 838 if(ic==1)then 839 d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 840 else 841 d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor) 842 end if 843 end do 844 end if 845 end if 846 end do 847 !$OMP END DO 848 end if 849 !$OMP END PARALLEL 850 end if 851 852 ! Compute <g|Vnl|c> (or derivatives) for each plane wave: 853 if (paw_opt/=3) then 854 !$OMP PARALLEL PRIVATE(ipw,ilmn,fdf,fdb,ffnl_dir1) 855 856 ! ------ 857 if (choice==1) then ! <g|Vnl|c> 858 !$OMP DO 859 do ipw=1,npw 860 ztab(ipw)=czero 861 do ilmn=1,nlmn 862 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 863 end do 864 end do 865 !$OMP END DO 866 867 ! ------ 868 else if (choice==2) then ! derivative w.r.t. atm. pos 869 !$OMP DO 870 do ipw=1,npw 871 ztab(ipw)=czero 872 do ilmn=1,nlmn 873 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp) 874 end do 875 ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw) 876 do ilmn=1,nlmn 877 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 878 end do 879 end do 880 !$OMP END DO 881 882 ! ------ 883 else if (choice==3) then ! derivative w.r.t. strain 884 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 885 if (idir<=3) then 886 !$OMP DO 887 do ipw=1,npw 888 ztab(ipw)=czero 889 do ilmn=1,nlmn 890 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) & 891 & *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp) & 892 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 893 end do 894 end do 895 !$OMP END DO 896 else 897 !$OMP DO 898 do ipw=1,npw 899 ztab(ipw)=czero 900 do ilmn=1,nlmn 901 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) & 902 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 903 end do 904 end do 905 !$OMP END DO 906 end if 907 908 909 ! ------ 910 else if (choice==5) then ! full derivative w.r.t. k 911 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 912 !$OMP DO 913 do ipw=1,npw 914 ztab(ipw)=czero 915 do ilmn=1,nlmn 916 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) & 917 & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 918 end do 919 end do 920 !$OMP END DO 921 922 ! ------ 923 else if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi> 924 !$OMP DO 925 do ipw=1,npw 926 ztab(ipw)=czero 927 do ilmn=1,nlmn 928 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 929 end do 930 end do 931 !$OMP END DO 932 933 ! ------ 934 else if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi> 935 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 936 !$OMP DO 937 do ipw=1,npw 938 ztab(ipw)=czero 939 do ilmn=1,nlmn 940 ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp) 941 end do 942 end do 943 !$OMP END DO 944 945 ! ------ 946 else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>V<dp/dk_(idir-1)|psi> 947 ! -<G|dp/dk_(idir-1)>V<dp/dk_(idir+1)|psi> 948 fdf = ffnl_dir_dat(2*idir-1) 949 fdb = ffnl_dir_dat(2*idir) 950 !$OMP DO 951 do ipw=1,npw 952 ztab(ipw)=czero 953 do ilmn=1,nlmn 954 ztab(ipw)=ztab(ipw) & 955 & +ffnl(ipw,fdf,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp) & 956 & -ffnl(ipw,fdb,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) 957 end do 958 end do 959 !$OMP END DO 960 961 ! ------ 962 else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k 963 !$OMP DO 964 do ipw=1,npw 965 ztab(ipw)=czero 966 do ilmn=1,nlmn 967 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp) 968 end do 969 ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw) 970 do ilmn=1,nlmn 971 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp) 972 end do 973 end do 974 !$OMP END DO 975 976 ! ------ 977 else if (choice==8) then ! full second order derivative w.r.t. k 978 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 979 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 980 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 981 !$OMP DO 982 do ipw=1,npw 983 ztab(ipw)=czero 984 do ilmn=1,nlmn 985 ztab(ipw)=ztab(ipw) & 986 & +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2, ilmn),kind=dp)& 987 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 988 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)& 989 & +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 990 end do 991 end do 992 !$OMP END DO 993 994 ! ------ 995 else if (choice==81) then 996 ! partial second order derivative w.r.t. k 997 ! full derivative w.r.t. k1, right derivative w.r.t. k2 998 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 999 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 1000 !$OMP DO 1001 do ipw=1,npw 1002 ztab(ipw)=czero 1003 do ilmn=1,nlmn 1004 ztab(ipw)=ztab(ipw) & 1005 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)& 1006 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp) 1007 end do 1008 end do 1009 !$OMP END DO 1010 1011 ! ------ 1012 else 1013 !$OMP WORKSHARE 1014 ztab(:)=czero 1015 !$OMP END WORKSHARE 1016 end if 1017 1018 1019 ! ------ 1020 !$OMP DO 1021 do ipw=1,npw 1022 ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp) 1023 vect(1,ipw+ipwshft)=vect(1,ipw+ipwshft)+real(ztab(ipw)) 1024 vect(2,ipw+ipwshft)=vect(2,ipw+ipwshft)+aimag(ztab(ipw)) 1025 end do 1026 !$OMP END DO 1027 1028 !$OMP END PARALLEL 1029 end if 1030 1031 ! Compute <g|S|c> (or derivatives) for each plane wave: 1032 if (paw_opt>=3) then 1033 !$OMP PARALLEL PRIVATE(ilmn,ipw,fdf,fdb,ffnl_dir1) 1034 1035 ! ------ 1036 if (choice==1) then ! <g|S|c> 1037 !$OMP DO 1038 do ipw=1,npw 1039 ztab(ipw)=czero 1040 do ilmn=1,nlmn 1041 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1042 end do 1043 end do 1044 !$OMP END DO 1045 1046 ! ------ 1047 else if (choice==2) then ! derivative w.r.t. atm. pos 1048 !$OMP DO 1049 do ipw=1,npw 1050 ztab(ipw)=czero 1051 do ilmn=1,nlmn 1052 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp) 1053 end do 1054 ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw) 1055 do ilmn=1,nlmn 1056 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 1057 end do 1058 end do 1059 !$OMP END DO 1060 1061 ! ------ 1062 else if (choice==3) then ! derivative w.r.t. strain 1063 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1064 if (idir<=3) then 1065 !$OMP DO 1066 do ipw=1,npw 1067 ztab(ipw)=czero 1068 do ilmn=1,nlmn 1069 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) & 1070 & *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)& 1071 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1072 end do 1073 end do 1074 !$OMP END DO 1075 else 1076 !$OMP DO 1077 do ipw=1,npw 1078 ztab(ipw)=czero 1079 do ilmn=1,nlmn 1080 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) & 1081 & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1082 end do 1083 end do 1084 !$OMP END DO 1085 end if 1086 1087 ! ------ 1088 else if (choice==5) then ! full derivative w.r.t. k 1089 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1090 !$OMP DO 1091 do ipw=1,npw 1092 ztab(ipw)=czero 1093 do ilmn=1,nlmn 1094 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) & 1095 & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1096 end do 1097 end do 1098 !$OMP END DO 1099 1100 ! ------ 1101 else if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi> 1102 !$OMP DO 1103 do ipw=1,npw 1104 ztab(ipw)=czero 1105 do ilmn=1,nlmn 1106 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 1107 end do 1108 end do 1109 !$OMP END DO 1110 1111 ! ------ 1112 else if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi> 1113 ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir 1114 !$OMP DO 1115 do ipw=1,npw 1116 ztab(ipw)=czero 1117 do ilmn=1,nlmn 1118 ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp) 1119 end do 1120 end do 1121 !$OMP END DO 1122 1123 ! ------ 1124 else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>S<dp/dk_(idir-1)|psi> - 1125 ! <G|dp/dk_(idir-1)>V<dp/dk_(idir+1)|psi> 1126 fdf = ffnl_dir_dat(2*idir-1) 1127 fdb = ffnl_dir_dat(2*idir) 1128 !$OMP DO 1129 do ipw=1,npw 1130 ztab(ipw)=czero 1131 do ilmn=1,nlmn 1132 ztab(ipw)=ztab(ipw) & 1133 & +ffnl(ipw,fdf,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp) & 1134 & -ffnl(ipw,fdb,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) 1135 end do 1136 end do 1137 !$OMP END DO 1138 1139 ! ------ 1140 else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k 1141 !$OMP DO 1142 do ipw=1,npw 1143 ztab(ipw)=czero 1144 do ilmn=1,nlmn 1145 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp) 1146 end do 1147 ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw) 1148 do ilmn=1,nlmn 1149 ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp) 1150 end do 1151 end do 1152 !$OMP END DO 1153 1154 ! ------ 1155 else if (choice==8) then ! full second order derivative w.r.t. k 1156 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1157 ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir) 1158 ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2)) 1159 !$OMP DO 1160 do ipw=1,npw 1161 ztab(ipw)=czero 1162 do ilmn=1,nlmn 1163 ztab(ipw)=ztab(ipw) & 1164 & +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx( gxfacs_(1, ilmn), gxfacs_(2, ilmn),kind=dp)& 1165 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 1166 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)& 1167 & +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 1168 end do 1169 end do 1170 !$OMP END DO 1171 1172 ! ------ 1173 else if (choice==81) then 1174 ! partial second order derivative w.r.t. k 1175 ! full derivative w.r.t. k1, right derivative w.r.t. k2 1176 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1177 ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir) 1178 !$OMP DO 1179 do ipw=1,npw 1180 ztab(ipw)=czero 1181 do ilmn=1,nlmn 1182 ztab(ipw)=ztab(ipw) & 1183 & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)& 1184 & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp) 1185 end do 1186 end do 1187 !$OMP END DO 1188 1189 ! ------ 1190 else 1191 !$OMP WORKSHARE 1192 ztab(:)=czero 1193 !$OMP END WORKSHARE 1194 end if 1195 1196 1197 ! ------ 1198 ! The OMP WORKSHARE directive doesn't have a good performance with Intel Compiler 1199 ! !$OMP WORKSHARE 1200 ! ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp) 1201 ! !$OMP END WORKSHARE 1202 ! !$OMP WORKSHARE 1203 ! svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:)) 1204 ! svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:)) 1205 ! !$OMP END WORKSHARE 1206 !$OMP DO 1207 do ipw=1,npw 1208 ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp) 1209 svect(1,ipw+ipwshft)=svect(1,ipw+ipwshft)+real(ztab(ipw)) 1210 svect(2,ipw+ipwshft)=svect(2,ipw+ipwshft)+aimag(ztab(ipw)) 1211 end do 1212 !$OMP END DO 1213 !$OMP END PARALLEL 1214 end if 1215 1216 ! End loop on atoms 1217 end do 1218 ! End loop on spinors 1219 end do 1220 1221 ! ========================================================================== 1222 end if 1223 1224 ABI_DEALLOCATE(ztab) 1225 1226 if (paw_opt/=3) then 1227 ABI_DEALLOCATE(gxfac_) 1228 if (choice>1) then 1229 ABI_DEALLOCATE(dgxdtfac_) 1230 end if 1231 if (choice==8) then 1232 ABI_DEALLOCATE(d2gxdtfac_) 1233 end if 1234 end if 1235 if (paw_opt>=3) then 1236 ABI_DEALLOCATE(gxfacs_) 1237 if (choice>1) then 1238 ABI_DEALLOCATE(dgxdtfacs_) 1239 end if 1240 if (choice==8) then 1241 ABI_DEALLOCATE(d2gxdtfacs_) 1242 end if 1243 end if 1244 1245 DBG_EXIT("COLL") 1246 1247 #if !defined HAVE_OPENMP 1248 !Fake use of unused variable 1249 if (.false.) write(std_out,*) ipw 1250 #endif 1251 1252 end subroutine opernlb_ylm