TABLE OF CONTENTS
ABINIT/opernl4b [ Functions ]
NAME
opernl4b
FUNCTION
Operate with the non-local part of the hamiltonian, from projected quantities to reciprocal space
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, 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
if(choice==2 .or choice==4 .or. choice==5) dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords or with respect to ddk if(choice==3) dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains ------ Taken away in beautification because unused MS ------- d2gxds2((2,mlang6,mincat,mproj) dummy argument here, not used ------------------------------------------------------------- ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere. gmet(3,3)=metric tensor for G vecs (in bohr**-2) gxa(2,mlang3,mincat,mproj)= projected scalars ia3=gives the number of the first atom in the subset presently treated idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1) indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln ispinor=1 or 2, gives the spinorial component of ffnl to be used ------ Taken away in beautification because unused MS ------- istwf_k=option parameter that describes the storage of wfs ------------------------------------------------------------- itypat = type of atom, needed for ffnl jproj(nlang)=number of projectors for each angular momentum kg_k(3,npw)=integer coords of planewaves in basis sphere kpg_k(npw,npkg)= (k+G) components and related data kpt(3)=real components of k point in terms of recip. translations lmnmax=max. number of (l,n) components over all type of psps matblk=dimension of the array ph3d mincat= maximum increment of atoms mlang3 = one of the dimensions of the array gxa mlang4 = dimension for dgxds ------ Taken away in beautification because unused MS ------- mlang6 = dimension for d2gxds2 ------------------------------------------------------------- mproj=maximum dimension for number of projection operators for each angular momentum for nonlocal pseudopotential ndgxdt=second dimension of dgxdt nincat = number of atoms in the subset here treated nkpg=second size of array kpg_k nffnl=third dimension of ffnl nlang = number of angular momenta to be treated = 1 + highest ang. mom. nloalg(3)=governs the choice of the algorithm for non-local operator. npw = number of plane waves in reciprocal space ntypat = number of type of atoms, dimension needed for ffnl gxa(2,mlang3,nincat,mproj)=modified projected scalars; NOTE that metric contractions have already been performed ph3d(2,npw,matblk)=three-dimensional phase factors
OUTPUT
vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.
NOTES
Operate with the non-local part of the hamiltonian for one type of atom, and within this given type of atom, for a subset of at most nincat atoms. This routine basically replaces getgla (gxa here is the former gla), except for the calculation of <G|dVnl/dk|C> or strain gradients. Present version decomposed according to iffkg opernl4a.f is from reciprocal space to projected quantities.
PARENTS
nonlop_pl
CHILDREN
dfpt_mkffkg
SOURCE
85 #if defined HAVE_CONFIG_H 86 #include "config.h" 87 #endif 88 89 #include "abi_common.h" 90 91 subroutine opernl4b(choice,dgxds,dgxdt,ffnl,gmet,gxa,& 92 & ia3,idir,indlmn,ispinor,itypat,jproj,kg_k,kpg_k,kpt,& 93 & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat,& 94 & nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) 95 96 use defs_basis 97 use m_profiling_abi 98 use m_errors 99 100 !This section has been created automatically by the script Abilint (TD). 101 !Do not modify the following lines by hand. 102 #undef ABI_FUNC 103 #define ABI_FUNC 'opernl4b' 104 use interfaces_66_nonlocal, except_this_one => opernl4b 105 !End of the abilint section 106 107 implicit none 108 109 !Arguments ------------------------------------ 110 !scalars 111 integer,intent(in) :: choice,ia3,idir,ispinor,itypat,lmnmax,matblk !,istwf_k 112 integer,intent(in) :: mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat !,mlang6 113 integer,intent(in) :: nkpg,nlang,npw,ntypat 114 !arrays 115 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 116 integer,intent(in) :: nloalg(3) 117 !real(dp),intent(in) :: d2gxds2(2,mlang6,mincat,mproj) 118 real(dp),intent(in) :: dgxds(2,mlang4,mincat,mproj) 119 real(dp),intent(in) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) 120 real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3) 121 real(dp),intent(in) :: gxa(2,mlang3,mincat,mproj),kpg_k(npw,nkpg),kpt(3) 122 real(dp),intent(in) :: ph3d(2,npw,matblk) 123 real(dp),intent(inout) :: vect(:,:) !vz_i 124 125 !Local variables------------------------------- 126 !scalars 127 integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,ig,ii,ilang,ilang2,ilang3 128 integer :: iproj,ipw,ipw1,ipw2,jffkg,jj,jump,mblkpw,nffkg 129 integer :: nffkgd,nffkge,nffkgk,nffkgs,nincpw,nproj,ntens,start 130 real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8 131 real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8 132 character(len=500) :: message 133 !arrays 134 integer,allocatable :: parity(:) 135 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 136 137 ! ************************************************************************* 138 139 !mblkpw sets the size of blocks of planewaves 140 mblkpw=NLO_MBLKPW 141 142 !jump governs, in fine, the use of registers in the most cpu 143 !time consuming part of the routine. Until now, jump=8 is the maximal value. 144 !The optimal value will be machine-dependent ! 145 jump=4 146 147 !Initialisation before blocking on the plane waves 148 149 !Set up dimension of kpgx and allocate 150 !ntens sets the maximum number of independent tensor components 151 !over all allowed angular momenta; need 20 for spdf for tensors 152 !up to rank 3; to handle stress tensor, need up to rank 5 153 ntens=1 154 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5)ntens=4 155 if(nlang>=3 .or. choice==3)ntens=10 156 if(nlang>=4 .or. (choice==3 .and. nlang>=2) )ntens=20 157 if(choice==3 .and. nlang>=3)ntens=35 158 if(choice==3 .and. nlang==4)ntens=56 159 160 !Set up second dimension of ffkg array, and allocate 161 nffkg=0; nffkge=0; nffkgd=0; nffkgk=0; nffkgs=0 162 do ilang=1,nlang 163 ! Get the number of projectors for that angular momentum 164 nproj=jproj(ilang) 165 ! If there is a non-local part, accumulate the number of vectors needed 166 if(nproj>0)then 167 ilang2=(ilang*(ilang+1))/2 168 nffkge=nffkge+nproj*ilang2 169 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 170 if(choice==2 .or. choice==4)nffkgd=nffkgd+ndgxdt*nproj*ilang2 171 if(choice==3)then 172 ilang3=((ilang+2)*(ilang+3))/2 173 nffkgs=nffkgs+nproj*ilang3 174 end if 175 end if 176 end do 177 nffkg=nffkge+nffkgd+nffkgs+nffkgk 178 179 !!$OMP PARALLEL DEFAULT(PRIVATE) & 180 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) & 181 !!$OMP SHARED(itypat,jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang3,mlang4,mproj) & 182 !!$OMP SHARED(ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) & 183 !!$OMP SHARED(jump,nffkgd,nffkgk,nffkgs,mblkpw,nffkg,nffkge,ntens) 184 185 ABI_ALLOCATE(ffkg,(nffkg,mblkpw)) 186 ABI_ALLOCATE(parity,(nffkg)) 187 ABI_ALLOCATE(kpgx,(mblkpw,ntens)) 188 ABI_ALLOCATE(scalars,(2,nffkg)) 189 ABI_ALLOCATE(teffv,(2,mblkpw)) 190 191 !Loop on subsets of plane waves (blocking) 192 !!$OMP DO 193 do ipw1=1,npw,mblkpw 194 195 ipw2=min(npw,ipw1+mblkpw-1) 196 nincpw=ipw2-ipw1+1 197 198 ! Initialize kpgx array related to tensors defined below 199 call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 200 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 201 & npw,ntens,ntypat,parity) 202 203 if (choice==1 .or. choice==2 .or. choice==3 .or. choice==5) then 204 ! Application of non-local part from projected scalars 205 ! back to reciprocal space ... 206 ! [this section merely computes terms which add to <G|Vnl|C>; 207 ! nothing here is needed when various gradients are being computed] 208 209 ! Loop on atoms 210 do ia=1,nincat 211 212 ! Compute the shift eventually needed to get the phases in ph3d 213 iaph3d=ia 214 if(nloalg(2)>0)iaph3d=ia+ia3-1 215 216 ! Transfer gxa (and eventually dgxdt) in scalars with different indexing 217 iffkg=0 218 iffkgk=nffkge*2 219 iffkgs=nffkge 220 do ilang=1,nlang 221 nproj=jproj(ilang) 222 if (nproj>0) then 223 ilang2=(ilang*(ilang+1))/2 224 ilang3=((ilang+2)*(ilang+3))/2 225 do iproj=1,nproj 226 do ii=1,ilang2 227 jj=ii+((ilang-1)*ilang*(ilang+1))/6 228 iffkg=iffkg+1 229 if(choice==1 .or. choice==3)then 230 scalars(1,iffkg)=gxa(1,jj,ia,iproj) 231 scalars(2,iffkg)=gxa(2,jj,ia,iproj) 232 else if (choice==2 .and. ndgxdt==1) then 233 jffkg=nffkge+iffkg 234 ! Pay attention to the use of reals and imaginary parts here ... 235 ! Also, the gxa and dgxdt arrays are switched, in order 236 ! to give the correct combination when multiplying ffkg, 237 ! see Eq.(53) of PRB55,10337(1997) 238 scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj) 239 scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj) 240 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 241 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 242 else if (choice==5) then 243 jffkg=nffkge+iffkg 244 ! The gxa and dgxdt arrays are switched, in order 245 ! to give the correct combination when multiplying ffkg, 246 scalars(1,jffkg)= gxa(1,jj,ia,iproj) 247 scalars(2,jffkg)= gxa(2,jj,ia,iproj) 248 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 249 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 250 end if 251 end do 252 if(choice==3) then 253 do ii=1,ilang3 254 iffkgs=iffkgs+1 255 jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 256 scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj) 257 scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj) 258 end do 259 end if 260 if(ilang>=2 .and. choice==5)then 261 do ii=1,((ilang-1)*ilang)/2 262 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 263 iffkgk=iffkgk+1 264 scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj) 265 scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj) 266 end do 267 end if 268 end do 269 end if 270 end do 271 272 ! DEBUG 273 ! write(std_out,*)' opernl4b, write scalars ' 274 ! do iffkg=1,nffkg 275 ! write(std_out,*)iffkg,scalars(1:2,iffkg) 276 ! end do 277 ! ENDDEBUG 278 279 ! ******* Entering the second critical part **************************** 280 281 ! First, treat small nffkg; send treat the loop needed for big nffkg; 282 ! finally treat the end of the loop needed for big nffkg 283 284 ! For the time being, the maximal jump allowed is 8. 285 286 ! 1) Here, treat small nffkg 287 if(nffkg<=jump)then 288 289 select case(nffkg) 290 291 case(1) 292 293 scr1=scalars(1,1) ; sci1=scalars(2,1) 294 ig=ipw1 295 do ipw=1,nincpw 296 ar=ffkg(1,ipw)*scr1 ; ai=ffkg(1,ipw)*sci1 297 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 298 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 299 ig=ig+1 300 end do 301 302 case(2) 303 304 scr1=scalars(1,1) ; sci1=scalars(2,1) 305 scr2=scalars(1,2) ; sci2=scalars(2,2) 306 ig=ipw1 307 do ipw=1,nincpw 308 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 309 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 310 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 311 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 312 ig=ig+1 313 end do 314 315 case(3) 316 317 scr1=scalars(1,1) ; sci1=scalars(2,1) 318 scr2=scalars(1,2) ; sci2=scalars(2,2) 319 scr3=scalars(1,3) ; sci3=scalars(2,3) 320 ig=ipw1 321 do ipw=1,nincpw 322 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 323 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 324 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 325 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 326 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 327 ig=ig+1 328 end do 329 330 case(4) 331 332 scr1=scalars(1,1) ; sci1=scalars(2,1) 333 scr2=scalars(1,2) ; sci2=scalars(2,2) 334 scr3=scalars(1,3) ; sci3=scalars(2,3) 335 scr4=scalars(1,4) ; sci4=scalars(2,4) 336 ig=ipw1 337 do ipw=1,nincpw 338 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 339 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 340 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 341 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 342 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 343 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 344 ig=ig+1 345 end do 346 347 case(5) 348 349 scr1=scalars(1,1) ; sci1=scalars(2,1) 350 scr2=scalars(1,2) ; sci2=scalars(2,2) 351 scr3=scalars(1,3) ; sci3=scalars(2,3) 352 scr4=scalars(1,4) ; sci4=scalars(2,4) 353 scr5=scalars(1,5) ; sci5=scalars(2,5) 354 ig=ipw1 355 do ipw=1,nincpw 356 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 357 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 358 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 359 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 360 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 361 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 362 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 363 ig=ig+1 364 end do 365 366 case(6) 367 368 scr1=scalars(1,1) ; sci1=scalars(2,1) 369 scr2=scalars(1,2) ; sci2=scalars(2,2) 370 scr3=scalars(1,3) ; sci3=scalars(2,3) 371 scr4=scalars(1,4) ; sci4=scalars(2,4) 372 scr5=scalars(1,5) ; sci5=scalars(2,5) 373 scr6=scalars(1,6) ; sci6=scalars(2,6) 374 ig=ipw1 375 do ipw=1,nincpw 376 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 377 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 378 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 379 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 380 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 381 ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6 382 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 383 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 384 ig=ig+1 385 end do 386 387 case(7) 388 389 scr1=scalars(1,1) ; sci1=scalars(2,1) 390 scr2=scalars(1,2) ; sci2=scalars(2,2) 391 scr3=scalars(1,3) ; sci3=scalars(2,3) 392 scr4=scalars(1,4) ; sci4=scalars(2,4) 393 scr5=scalars(1,5) ; sci5=scalars(2,5) 394 scr6=scalars(1,6) ; sci6=scalars(2,6) 395 scr7=scalars(1,7) ; sci7=scalars(2,7) 396 ig=ipw1 397 do ipw=1,nincpw 398 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 399 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 400 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 401 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 402 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 403 ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6 404 ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7 405 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 406 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 407 ig=ig+1 408 end do 409 410 case(8) 411 412 scr1=scalars(1,1) ; sci1=scalars(2,1) 413 scr2=scalars(1,2) ; sci2=scalars(2,2) 414 scr3=scalars(1,3) ; sci3=scalars(2,3) 415 scr4=scalars(1,4) ; sci4=scalars(2,4) 416 scr5=scalars(1,5) ; sci5=scalars(2,5) 417 scr6=scalars(1,6) ; sci6=scalars(2,6) 418 scr7=scalars(1,7) ; sci7=scalars(2,7) 419 scr8=scalars(1,8) ; sci8=scalars(2,8) 420 ig=ipw1 421 do ipw=1,nincpw 422 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 423 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 424 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 425 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 426 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 427 ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6 428 ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7 429 ar=ar+ffkg(8,ipw)*scr8 ; ai=ai+ffkg(8,ipw)*sci8 430 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 431 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 432 ig=ig+1 433 end do 434 435 436 end select 437 438 else 439 440 do ipw=1,nincpw 441 teffv(1,ipw)=0.0d0 ; teffv(2,ipw)=0.0d0 442 end do 443 444 ! 2) Here treart the loop for big nffkg 445 do start=1,nffkg-jump,jump 446 chunk=min(jump,nffkg-jump-start+1) 447 448 select case(chunk) 449 450 case(1) 451 452 scr1=scalars(1,start) ; sci1=scalars(2,start) 453 do ipw=1,nincpw 454 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 455 ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1 456 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 457 end do 458 459 case(2) 460 461 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 462 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 463 do ipw=1,nincpw 464 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 465 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 466 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 467 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 468 end do 469 470 case(3) 471 472 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 473 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 474 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 475 do ipw=1,nincpw 476 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 477 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 478 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 479 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 480 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 481 end do 482 483 case(4) 484 485 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 486 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 487 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 488 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 489 do ipw=1,nincpw 490 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 491 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 492 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 493 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 494 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 495 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 496 end do 497 498 case(5) 499 500 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 501 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 502 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 503 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 504 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 505 do ipw=1,nincpw 506 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 507 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 508 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 509 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 510 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 511 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 512 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 513 end do 514 515 case(6) 516 517 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 518 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 519 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 520 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 521 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 522 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 523 do ipw=1,nincpw 524 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 525 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 526 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 527 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 528 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 529 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 530 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 531 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 532 end do 533 534 case(7) 535 536 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 537 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 538 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 539 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 540 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 541 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 542 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 543 do ipw=1,nincpw 544 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 545 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 546 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 547 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 548 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 549 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 550 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 551 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 552 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 553 end do 554 555 case(8) 556 557 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 558 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 559 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 560 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 561 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 562 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 563 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 564 scr8=scalars(1,start+7) ; sci8=scalars(2,start+7) 565 do ipw=1,nincpw 566 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 567 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 568 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 569 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 570 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 571 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 572 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 573 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 574 ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8 575 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 576 end do 577 578 end select 579 580 end do 581 582 ! 3) Treat the end of the loops 583 584 start=nffkg-jump+1 585 586 select case(jump) 587 588 case(1) 589 590 scr1=scalars(1,start) ; sci1=scalars(2,start) 591 ig=ipw1 592 do ipw=1,nincpw 593 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 594 ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1 595 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 596 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 597 ig=ig+1 598 end do 599 600 case(2) 601 602 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 603 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 604 ig=ipw1 605 do ipw=1,nincpw 606 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 607 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 608 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 609 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 610 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 611 ig=ig+1 612 end do 613 614 case(3) 615 616 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 617 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 618 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 619 ig=ipw1 620 do ipw=1,nincpw 621 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 622 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 623 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 624 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 625 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 626 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 627 ig=ig+1 628 end do 629 630 case(4) 631 632 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 633 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 634 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 635 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 636 ig=ipw1 637 do ipw=1,nincpw 638 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 639 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 640 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 641 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 642 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 643 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 644 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 645 ig=ig+1 646 end do 647 648 case(5) 649 650 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 651 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 652 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 653 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 654 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 655 ig=ipw1 656 do ipw=1,nincpw 657 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 658 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 659 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 660 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 661 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 662 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 663 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 664 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 665 ig=ig+1 666 end do 667 668 case(6) 669 670 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 671 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 672 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 673 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 674 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 675 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 676 ig=ipw1 677 do ipw=1,nincpw 678 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 679 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 680 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 681 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 682 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 683 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 684 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 685 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 686 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 687 ig=ig+1 688 end do 689 690 case(7) 691 692 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 693 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 694 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 695 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 696 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 697 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 698 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 699 ig=ipw1 700 do ipw=1,nincpw 701 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 702 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 703 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 704 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 705 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 706 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 707 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 708 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 709 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 710 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 711 ig=ig+1 712 end do 713 714 case(8) 715 716 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 717 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 718 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 719 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 720 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 721 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 722 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 723 scr8=scalars(1,start+7) ; sci8=scalars(2,start+7) 724 ig=ipw1 725 do ipw=1,nincpw 726 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 727 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 728 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 729 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 730 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 731 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 732 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 733 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 734 ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8 735 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 736 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 737 ig=ig+1 738 end do 739 740 end select 741 742 ! End if statement for small or big nffkg 743 end if 744 745 ! ******* Leaving the critical part ********************************* 746 747 ! End loop on atoms 748 end do 749 750 ! End choice==1 or choice==2 or choice==3 751 else 752 ! Problem: choice does not make sense 753 write(message,'(a,i0,a)' )' Input choice=',choice,' not allowed. ' 754 MSG_BUG(message) 755 end if 756 757 ! End loop on blocks of planewaves 758 end do 759 !!$OMP END DO 760 761 ABI_DEALLOCATE(ffkg) 762 ABI_DEALLOCATE(kpgx) 763 ABI_DEALLOCATE(parity) 764 ABI_DEALLOCATE(scalars) 765 ABI_DEALLOCATE(teffv) 766 !!$OMP END PARALLEL 767 768 end subroutine opernl4b