TABLE OF CONTENTS
ABINIT/m_opernl [ Modules ]
NAME
m_opernl
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 ABINIT group (DCA, XG, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_opernl 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 28 use m_mkffkg, only : mkffkg, dfpt_mkffkg 29 30 implicit none 31 32 private
ABINIT/opernl2 [ Functions ]
NAME
opernl2
FUNCTION
Operate with the non-local part of the hamiltonian, either from reciprocal space to projected quantities (sign=1), or from projected quantities to reciprocal space (sign=-1)
INPUTS
if(sign==-1 .and. (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(sign==-1 .and. choice==3) dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere. gmet(3,3)=metric tensor for G vecs (in bohr**-2) 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 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 mlang1 = dimensions for dgxdis1 mlang3 = one of the dimensions of the array gxa mlang4 = dimension for dgxds mlang5 = dimensions for dgxdis2 mlang6 = dimension for d2gxds2 mproj=maximum dimension for number of projection operators for each angular momentum for nonlocal pseudopotential ndgxdt=second dimension of dgxdt nffnl=second dimension of ffnl nincat = number of atoms in the subset here treated nkpg=second size of array kpg_k 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 sign : if 1, go from reciprocal space to projected scalars, if -1, go from projected scalars to reciprocal space. if(sign==1), vect(2*npw)=starting vector in reciprocal space if(sign==-1) gxa(2,mlang3,nincat,mproj)=modified projected scalars; NOTE that metric contractions have already been performed on the arrays gxa if sign=-1 ph3d(2,npw,matblk)=three-dimensional phase factors
OUTPUT
if(sign==1) gxa(2,mlang3,mincat,mproj)= projected scalars if(sign==1 .and. (choice==2 .or choice==4 .or. choice==5 .or. choice==23)) dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords or with respect to ddk if(sign==1 .and. (choice==3 .or. choice==23)) dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains if(sign==1 .and. choice==6) dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars wrt coord. indexed for internal strain d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars wrt strain and coord d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars wrt strains if(sign==-1) 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.
SOURCE
125 subroutine opernl2(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 126 & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,& 127 & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 128 & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,& 129 & ntypat,ph3d,sign,vect) 130 131 132 !Arguments ------------------------------------ 133 !scalars 134 integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk 135 integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt 136 integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign 137 !arrays 138 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 139 integer,intent(in) :: nloalg(3) 140 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 141 real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk) 142 real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) 143 real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) 144 real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:) 145 real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj) 146 real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj) 147 real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj) 148 149 !Local variables------------------------------- 150 !scalars 151 integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3 152 integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw 153 integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj,ntens 154 real(dp),parameter :: two_pi2=two_pi**2 155 character(len=500) :: message 156 !arrays 157 integer,allocatable :: parity(:) 158 !real(dp) :: tsec(2) 159 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 160 161 ! ************************************************************************* 162 163 !call wrtout(std_out,"in opernl2","COLL") 164 165 !mblkpw sets the size of blocks of planewaves 166 mblkpw=NLO_MBLKPW 167 168 !Get the actual maximum number of projectors 169 mmproj=maxval(indlmn(3,:,itypat)) 170 171 !Initialisation before blocking on the plane waves 172 173 if (sign==1) then 174 ! Put projected scalars to zero 175 gxa(:,:,:,1:mmproj)=0.0d0 176 if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0 177 if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0 178 if (choice==6) then 179 dgxdis(:,:,:,1:mmproj)=0.0d0 180 d2gxdis(:,:,:,1:mmproj)=0.0d0 181 d2gxds2(:,:,:,1:mmproj)=0.0d0 182 end if 183 end if 184 185 !Set up dimension of kpgx and allocate 186 !ntens sets the maximum number of independent tensor components 187 !over all allowed angular momenta; need 20 for spdf for tensors 188 !up to rank 3; to handle stress tensor, need up to rank 5 189 ntens=1 190 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4 191 if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10 192 if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20 193 if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35 194 if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56 195 if(choice==6 .and. nlang>=3)ntens=84 196 if(choice==6 .and. nlang==4)ntens=120 197 198 !Set up second dimension of ffkg array, and allocate 199 nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0 200 do ilang=1,nlang 201 ! Get the number of projectors for that angular momentum 202 nproj=jproj(ilang) 203 ! If there is a non-local part, accumulate the number of vectors needed 204 ! The variables ilang below are the number of independent tensors of 205 ! various ranks, the variable names being more historical than logical. 206 ! ilang2=number of rank ilang-1 207 ! ilang3=number of rank ilang+1 208 ! ilang4=number of rank ilang 209 ! ilang5=number of rank ilang+2 210 ! ilang6=number of rank ilang+3 211 if(nproj>0)then 212 ilang2=(ilang*(ilang+1))/2 213 nffkge=nffkge+nproj*ilang2 214 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 215 if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2 216 if(choice==3 .or. choice==6 .or. choice==23)then 217 ilang3=((ilang+2)*(ilang+3))/2 218 nffkgs=nffkgs+nproj*ilang3 219 end if 220 if(choice==6)then 221 ilang4=((ilang+1)*(ilang+2))/2 222 ilang5=((ilang+3)*(ilang+4))/2 223 ilang6=((ilang+4)*(ilang+5))/2 224 nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6) 225 end if 226 end if 227 end do 228 nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk 229 230 !Loop on subsets of plane waves (blocking) 231 !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a sigfault with nthreads==1 232 !Feb 16 2012: The code does not crash anymore but it's not efficient. 233 ! 234 !!$OMP PARALLEL DEFAULT(PRIVATE) & 235 !!$OMP SHARED(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) & 236 !!$OMP SHARED(ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) & 237 !!$OMP SHARED(jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) & 238 !!$OMP SHARED(mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) & 239 !!$OMP SHARED(ntypat,ph3d,sign,vect) & 240 !!$OMP SHARED(mblkpw,nffkg,nffkgd,nffkge,nffkgs,ntens,mmproj) 241 242 ABI_MALLOC(ffkg,(mblkpw,nffkg)) 243 ABI_MALLOC(parity,(nffkg)) 244 ABI_MALLOC(kpgx,(mblkpw,ntens)) 245 ABI_MALLOC(scalars,(2,nffkg)) 246 ABI_MALLOC(teffv,(2,mblkpw)) 247 !!$OMP DO 248 do ipw1=1,npw,mblkpw 249 250 ipw2=min(npw,ipw1+mblkpw-1) 251 nincpw=ipw2-ipw1+1 252 253 ! call timab(74+choice,1,tsec) 254 255 ! Initialize kpgx array related to tensors defined below 256 call mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,kg_k,& 257 & kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 258 & npw,ntens,ntypat,parity) 259 260 ! call timab(74+choice,2,tsec) 261 262 ! Now treat the different signs 263 if (sign==1) then 264 265 do ia=1,nincat 266 267 ! Compute the shift eventually needed to get the phases in ph3d 268 iaph3d=ia 269 if(nloalg(2)>0)iaph3d=ia+ia3-1 270 271 ! ******* Entering the first time-consuming part of the routine ******* 272 273 ! Multiply by the phase factor 274 ! This allows to be left with only real operations, 275 ! that are performed in the most inner loops 276 ig=ipw1 277 do ipw=1,nincpw 278 teffv(1,ipw)=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 279 teffv(2,ipw)=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 280 ig=ig+1 281 end do 282 283 do iffkg=1,nffkg 284 scalars(1,iffkg)=0.0d0 285 scalars(2,iffkg)=0.0d0 286 do ipw=1,nincpw 287 scalars(1,iffkg)=scalars(1,iffkg)+teffv(1,ipw)*ffkg(ipw,iffkg) 288 scalars(2,iffkg)=scalars(2,iffkg)+teffv(2,ipw)*ffkg(ipw,iffkg) 289 end do 290 end do 291 292 ! ******* Leaving the critical part ********************************* 293 294 ! DEBUG 295 ! write(std_out,*)' opernl2 : scalars' 296 ! do iffkg=1,10 297 ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg) 298 ! end do 299 ! stop 300 ! ENDDEBUG 301 302 if(istwf_k>=2)then 303 ! Impose parity of resulting scalar (this operation could be 304 ! replaced by direct saving of CPU time in the preceeding section) 305 do iffkg=1,nffkg 306 scalars(parity(iffkg),iffkg)=0.0d0 307 end do 308 end if 309 310 iffkg=0 311 iffkgs=nffkge+nffkgd 312 iffkgs2=nffkge+nffkgs 313 iffkgk=nffkge*2 314 do ilang=1,nlang 315 nproj=jproj(ilang) 316 if(nproj>0)then 317 ! ilang2 is the number of independent tensor components 318 ! for symmetric tensor of rank ilang-1 319 ilang2=(ilang*(ilang+1))/2 320 321 ! Loop over projectors 322 do iproj=1,nproj 323 ! Multiply by the k+G factors (tensors of various rank) 324 do ii=1,ilang2 325 ! Get the starting address for the relevant tensor 326 jj=ii+((ilang-1)*ilang*(ilang+1))/6 327 iffkg=iffkg+1 328 ! !$OMP CRITICAL (OPERNL2_1) 329 gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg) 330 gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg) 331 ! !$OMP END CRITICAL (OPERNL2_1) 332 ! Now, compute gradients, if needed. 333 if ((choice==2.or.choice==23) .and. ndgxdt==3) then 334 do mu=1,3 335 jffkg=nffkge+(iffkg-1)*3+mu 336 ! Pay attention to the use of reals and imaginary parts here ... 337 ! !$OMP CRITICAL (OPERNL2_2) 338 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 339 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 340 ! !$OMP END CRITICAL (OPERNL2_2) 341 end do 342 end if 343 if (choice==2 .and. ndgxdt==1) then 344 jffkg=nffkge+iffkg 345 ! Pay attention to the use of reals and imaginary parts here ... 346 ! !$OMP CRITICAL (OPERNL2_3) 347 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg) 348 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg) 349 ! !$OMP END CRITICAL (OPERNL2_3) 350 end if 351 if (choice==4) then 352 do mu=1,3 353 jffkg=nffkge+(iffkg-1)*9+mu 354 ! Pay attention to the use of reals and imaginary parts here ... 355 ! !$OMP CRITICAL (OPERNL2_4) 356 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 357 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 358 ! !$OMP END CRITICAL (OPERNL2_4) 359 end do 360 do mu=4,9 361 jffkg=nffkge+(iffkg-1)*9+mu 362 ! Pay attention to the use of reals and imaginary parts here ... 363 ! Also, note the multiplication by (2 pi)**2 364 ! !$OMP CRITICAL (OPERNL2_5) 365 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg) 366 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg) 367 ! !$OMP END CRITICAL (OPERNL2_5) 368 end do 369 end if 370 ! End loop on ii=1,ilang2 371 end do 372 373 if (choice==3 .or. choice==6 .or. choice==23) then 374 ! Compute additional tensors related to strain gradients 375 ! ilang3 is number of unique tensor components of rank ilang+1 376 ilang3=((ilang+2)*(ilang+3))/2 377 jjs=((ilang+1)*(ilang+2)*(ilang+3))/6 378 ! Compute strain gradient tensor components 379 do ii=1,ilang3 380 ! Note that iffkgs is also used by ddk and 2nd derivative parts 381 iffkgs=iffkgs+1 382 jj=ii+jjs 383 ! !$OMP CRITICAL (OPERNL2_6) 384 dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs) 385 dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs) 386 ! !$OMP END CRITICAL (OPERNL2_6) 387 end do 388 end if 389 390 if (choice==6) then 391 ! Compute additional tensors related to strain 2nd derivatives 392 ! and internal strain derivatives 393 ! ilang6 is number of unique tensor components of rank ilang+3 394 ilang6=((ilang+4)*(ilang+5))/2 395 jjs=((ilang+3)*(ilang+4)*(ilang+5))/6 396 ! Compute strain gradient tensor components 397 do ii=1,ilang6 398 iffkgs2=iffkgs2+1 399 jj=ii+jjs 400 ! !$OMP CRITICAL (OPERNL2_6) 401 d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2) 402 d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2) 403 ! !$OMP END CRITICAL (OPERNL2_6) 404 end do 405 406 ! ilang4 is number of unique tensor components of rank ilang 407 ilang4=((ilang+1)*(ilang+2))/2 408 jjs=((ilang)*(ilang+1)*(ilang+2))/6 409 ! Compute internal strain gradient tensor components 410 do ii=1,ilang4 411 iffkgs2=iffkgs2+1 412 jj=ii+jjs 413 ! !$OMP CRITICAL (OPERNL2_6) 414 ! Pay attention to the use of reals and imaginary parts here ... 415 dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2) 416 dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2) 417 ! !$OMP END CRITICAL (OPERNL2_6) 418 end do 419 420 ! ilang5 is number of unique tensor components of rank ilang+2 421 ilang5=((ilang+3)*(ilang+4))/2 422 jjs=((ilang+2)*(ilang+3)*(ilang+4))/6 423 ! Compute internal strain gradient tensor components 424 do ii=1,ilang5 425 iffkgs2=iffkgs2+1 426 jj=ii+jjs 427 ! !$OMP CRITICAL (OPERNL2_6) 428 ! Pay attention to the use of reals and imaginary parts here ... 429 d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2) 430 d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2) 431 ! !$OMP END CRITICAL (OPERNL2_6) 432 end do 433 end if ! choice==6 434 435 if (choice==5) then 436 ! Compute additional tensors related to ddk with ffnl(:,2,...) 437 ilangx=(ilang*(ilang+1))/2 438 jjs=((ilang-1)*ilang*(ilang+1))/6 439 do ii=1,ilangx 440 ! Note that iffkgs is also used by strain part 441 iffkgs=iffkgs+1 442 jj=ii+jjs 443 ! !$OMP CRITICAL (OPERNL2_7) 444 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs) 445 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs) 446 ! !$OMP END CRITICAL (OPERNL2_7) 447 end do 448 ! Compute additional tensors related to ddk with ffnl(:,1,...) 449 if(ilang>=2)then 450 ilangx=((ilang-1)*ilang)/2 451 jjs=((ilang-2)*(ilang-1)*ilang)/6 452 do ii=1,ilangx 453 iffkgk=iffkgk+1 454 jj=ii+jjs 455 ! !$OMP CRITICAL (OPERNL2_8) 456 dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk) 457 dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk) 458 ! !$OMP END CRITICAL (OPERNL2_8) 459 end do 460 end if 461 end if 462 463 ! End projector loop 464 end do 465 466 ! End condition of non-zero projectors 467 end if 468 469 ! End angular momentum loop 470 end do 471 472 ! End loop on atoms 473 end do 474 475 else if (sign==-1) then 476 ! Application of non-local part from projected scalars 477 ! back to reciprocal space ... 478 ! [this section merely computes terms which add to <G|Vnl|C>; 479 ! nothing here is needed when various gradients are being computed] 480 481 ! Loop on atoms 482 do ia=1,nincat 483 484 ! Compute the shift eventually needed to get the phases in ph3d 485 iaph3d=ia 486 if(nloalg(2)>0)iaph3d=ia+ia3-1 487 488 ! Transfer gxa (and eventually dgxdt) in scalars with different indexing 489 iffkg=0 490 iffkgk=nffkge*2 491 iffkgs=nffkge 492 do ilang=1,nlang 493 nproj=jproj(ilang) 494 if (nproj>0) then 495 ilang2=(ilang*(ilang+1))/2 496 ilang3=((ilang+2)*(ilang+3))/2 497 do iproj=1,nproj 498 do ii=1,ilang2 499 jj=ii+((ilang-1)*ilang*(ilang+1))/6 500 iffkg=iffkg+1 501 if(choice==1 .or. choice==3)then 502 scalars(1,iffkg)=gxa(1,jj,ia,iproj) 503 scalars(2,iffkg)=gxa(2,jj,ia,iproj) 504 else if (choice==2 .and. ndgxdt==1) then 505 jffkg=nffkge+iffkg 506 ! Pay attention to the use of reals and imaginary parts here ... 507 ! Also, the gxa and dgxdt arrays are switched, in order 508 ! to give the correct combination when multiplying ffkg, 509 ! see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]] 510 scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj) 511 scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj) 512 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 513 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 514 else if (choice==5) then 515 jffkg=nffkge+iffkg 516 ! The gxa and dgxdt arrays are switched, in order 517 ! to give the correct combination when multiplying ffkg, 518 scalars(1,jffkg)= gxa(1,jj,ia,iproj) 519 scalars(2,jffkg)= gxa(2,jj,ia,iproj) 520 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 521 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 522 end if 523 end do 524 if(choice==3) then 525 do ii=1,ilang3 526 iffkgs=iffkgs+1 527 jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 528 scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj) 529 scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj) 530 end do 531 end if 532 if(ilang>=2 .and. choice==5)then 533 do ii=1,((ilang-1)*ilang)/2 534 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 535 iffkgk=iffkgk+1 536 scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj) 537 scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj) 538 end do 539 end if 540 end do 541 end if 542 end do 543 544 ! DEBUG 545 ! if(choice==5)then 546 ! write(std_out,*)' opernl2 : write gxa(:,...) array ' 547 ! do ii=1,10 548 ! write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1) 549 ! end do 550 ! write(std_out,*)' opernl2 : write dgxdt(:,1,...) array ' 551 ! do ii=1,10 552 ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1) 553 ! end do 554 ! write(std_out,*)' opernl2 : write dgxdt(:,2,...) array ' 555 ! do ii=1,4 556 ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1) 557 ! end do 558 ! end if 559 560 ! do iffkg=1,nffkg 561 ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg) 562 ! end do 563 ! stop 564 ! ENDDEBUG 565 566 ! ******* Entering the critical part ********************************* 567 568 do ipw=1,nincpw 569 teffv(1,ipw)=0.0d0 570 teffv(2,ipw)=0.0d0 571 end do 572 do iffkg=1,nffkg 573 do ipw=1,nincpw 574 teffv(1,ipw)=teffv(1,ipw)+ffkg(ipw,iffkg)*scalars(1,iffkg) 575 teffv(2,ipw)=teffv(2,ipw)+ffkg(ipw,iffkg)*scalars(2,iffkg) 576 end do 577 end do 578 ! Multiplication by the complex conjugate of the phase 579 ig=ipw1 580 do ipw=1,nincpw 581 vect(1,ig)=vect(1,ig)+teffv(1,ipw)*ph3d(1,ig,iaph3d)+teffv(2,ipw)*ph3d(2,ig,iaph3d) 582 vect(2,ig)=vect(2,ig)+teffv(2,ipw)*ph3d(1,ig,iaph3d)-teffv(1,ipw)*ph3d(2,ig,iaph3d) 583 ig=ig+1 584 end do 585 586 ! ******* Leaving the critical part ********************************* 587 588 ! End loop on atoms 589 end do 590 591 ! End sign==-1 592 else 593 594 ! Problem: sign and/or choice do not make sense 595 write(message,'(a,2i10,a,a)')& 596 & ' Input sign, choice=',sign,choice,ch10,& 597 & ' Are not compatible or allowed. ' 598 ABI_BUG(message) 599 end if 600 601 ! End loop on blocks of planewaves 602 end do 603 !!$OMP END DO 604 ABI_FREE(ffkg) 605 ABI_FREE(kpgx) 606 ABI_FREE(parity) 607 ABI_FREE(scalars) 608 ABI_FREE(teffv) 609 !!$OMP END PARALLEL 610 611 612 !DEBUG 613 !if(choice==5)then 614 !write(std_out,*)'opernl2 : write vect(2*npw)' 615 !do ipw=1,2 616 !write(std_out,*)ipw,vect(1:2,ipw) 617 !end do 618 !write(std_out,*)'opernl2 : write ph3d' 619 !do ipw=1,npw 620 !write(std_out,*)ipw,ph3d(1:2,ipw,1) 621 !end do 622 !write(std_out,*)' opernl2 : write gxa array ' 623 !write(std_out,*)' ang mom ,ia ' 624 !do iproj=1,mproj 625 !do ia=1,1 626 !do ii=1,3 627 !do ii=1,10 628 !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1) 629 !end do 630 !end do 631 !end do 632 !end if 633 !if(choice==5)then 634 !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array ' 635 !do ii=1,10 636 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1) 637 !end do 638 !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array ' 639 !do ii=1,4 640 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1) 641 !end do 642 !stop 643 !end if 644 !ENDDEBUG 645 646 end subroutine opernl2
ABINIT/opernl3 [ Functions ]
NAME
opernl3
FUNCTION
Operate with the non-local part of the hamiltonian, either from reciprocal space to projected quantities (sign=1), or from projected quantities to reciprocal space (sign=-1)
INPUTS
if(sign==-1 .and. (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(sign==-1 .and. choice==3) dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere. gmet(3,3)=metric tensor for G vecs (in bohr**-2) 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 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 mlang1 = dimensions for dgxdis1 mlang3 = one of the dimensions of the array gxa mlang4 = dimension for dgxds mlang5 = dimensions for dgxdis2 mlang6 = dimension for d2gxds2 mproj=maximum dimension for number of projection operators for each angular momentum for nonlocal pseudopotential ndgxdt=second dimension of dgxdt nffnl=third dimension of ffnl nincat = number of atoms in the subset here treated nkpg=second size of array kpg_k 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 sign : if 1, go from reciprocal space to projected scalars, if -1, go from projected scalars to reciprocal space. if(sign==1), vect(2*npw)=starting vector in reciprocal space if(sign==-1) gxa(2,mlang3,nincat,mproj)=modified projected scalars; NOTE that metric contractions have already been performed on the arrays gxa if sign=-1 ph3d(2,npw,matblk)=three-dimensional phase factors
OUTPUT
if(sign==1) gxa(2,mlang3,mincat,mproj)= projected scalars if(sign==1 .and. (choice==2 .or choice==4 .or. choice==5 .or. choice==23)) dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords or with respect to ddk if(sign==1 .and. (choice==3 .or. choice==23)) dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains if(sign==1 .and. choice==6) dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars wrt coord. indexed for internal strain d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars wrt strain and coord d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars wrt strains if(sign==-1) 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.
SOURCE
732 subroutine opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 733 & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,& 734 & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 735 & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,& 736 & ntypat,ph3d,sign,vect) 737 738 !Arguments ------------------------------------ 739 !scalars 740 integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk 741 integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt 742 integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign 743 !arrays 744 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 745 integer,intent(in) :: nloalg(3) 746 real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3) 747 real(dp),intent(in) :: kpg_k(npw,nkpg),kpt(3),ph3d(2,npw,matblk) 748 real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) 749 real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) 750 real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:) 751 real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj) 752 real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj) 753 real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj) 754 755 !Local variables------------------------------- 756 !ntens sets the maximum number of independent tensor components 757 !over all allowed angular momenta; need 20 for spdf for tensors 758 !up to rank 3; to handle stress tensor, need up to rank 5 759 !to handle strain 2DER, need up to rank 7 760 !scalars 761 integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3 762 integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw 763 integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj 764 integer :: ntens 765 real(dp) :: ai,ar 766 real(dp),parameter :: two_pi2=two_pi**2 767 character(len=500) :: message 768 !arrays 769 integer,allocatable :: parity(:) 770 ! real(dp) :: tsec(2) 771 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 772 773 ! ************************************************************************* 774 775 !call wrtout(std_out,"in opernl3","COLL") 776 777 !mblkpw sets the size of blocks of planewaves 778 mblkpw=NLO_MBLKPW 779 780 !Get the actual maximum number of projectors 781 mmproj=maxval(indlmn(3,:,itypat)) 782 783 !Initialisation before blocking on the plane waves 784 785 if (sign==1) then 786 ! Put projected scalars to zero 787 gxa(:,:,:,1:mmproj)=0.0d0 788 if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0 789 if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0 790 if (choice==6) then 791 dgxdis(:,:,:,1:mmproj)=0.0d0 792 d2gxdis(:,:,:,1:mmproj)=0.0d0 793 d2gxds2(:,:,:,1:mmproj)=0.0d0 794 end if 795 end if 796 797 !Set up dimension of kpgx and allocate 798 !ntens sets the maximum number of independent tensor components 799 !over all allowed angular momenta; need 20 for spdf for tensors 800 !up to rank 3; to handle stress tensor, need up to rank 5 801 ntens=1 802 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4 803 if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10 804 if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20 805 if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35 806 if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56 807 if(choice==6 .and. nlang>=3)ntens=84 808 if(choice==6 .and. nlang==4)ntens=120 809 810 !Set up second dimension of ffkg array, and allocate 811 nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0 812 do ilang=1,nlang 813 ! Get the number of projectors for that angular momentum 814 nproj=jproj(ilang) 815 ! If there is a non-local part, accumulate the number of vectors needed 816 ! The variables ilang below are the number of independent tensors of 817 ! various ranks, the variable names being more historical than logical. 818 ! ilang2=number of rank ilang-1 819 ! ilang3=number of rank ilang+1 820 ! ilang4=number of rank ilang 821 ! ilang5=number of rank ilang+2 822 ! ilang6=number of rank ilang+3 823 if(nproj>0)then 824 ilang2=(ilang*(ilang+1))/2 825 nffkge=nffkge+nproj*ilang2 826 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 827 if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2 828 if(choice==3 .or. choice==6 .or. choice==23)then 829 ilang3=((ilang+2)*(ilang+3))/2 830 nffkgs=nffkgs+nproj*ilang3 831 end if 832 if(choice==6)then 833 ilang4=((ilang+1)*(ilang+2))/2 834 ilang5=((ilang+3)*(ilang+4))/2 835 ilang6=((ilang+4)*(ilang+5))/2 836 nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6) 837 end if 838 end if 839 end do 840 nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk 841 842 !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a 843 !sigfault with nthreads==1 844 ! 845 !Loop on subsets of plane waves (blocking) 846 !!$OMP PARALLEL DEFAULT(PRIVATE) & 847 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) & 848 !!$OMP SHARED(istwf_k,itypat,jproj,kg_k,kpg_k,kpt,lmnmax,mblkpw,mproj) & 849 !!$OMP SHARED(ndgxdt,nffkg,nffkgd,nffkge,nffkgs,nincat,nkpg,nlang) & 850 !!$OMP SHARED(nloalg,ph3d,npw,ntens,ntypat,sign,vect) 851 852 ABI_MALLOC(ffkg,(nffkg,mblkpw)) 853 ABI_MALLOC(parity,(nffkg)) 854 ABI_MALLOC(kpgx,(mblkpw,ntens)) 855 ABI_MALLOC(scalars,(2,nffkg)) 856 ABI_MALLOC(teffv,(2,mblkpw)) 857 !!$OMP DO 858 do ipw1=1,npw,mblkpw 859 860 ipw2=min(npw,ipw1+mblkpw-1) 861 nincpw=ipw2-ipw1+1 862 863 ! Initialize kpgx array related to tensors defined below 864 call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 865 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 866 & npw,ntens,ntypat,parity) 867 868 ! call timab(74+choice,2,tsec) 869 870 ! Now treat the different signs 871 if (sign==1) then 872 873 do ia=1,nincat 874 875 ! Compute the shift eventually needed to get the phases in ph3d 876 iaph3d=ia 877 if(nloalg(2)>0)iaph3d=ia+ia3-1 878 879 do iffkg=1,nffkg 880 scalars(1,iffkg)=0.0d0 881 scalars(2,iffkg)=0.0d0 882 end do 883 884 ! ******* Entering the first time-consuming part of the routine ******* 885 886 ! Note : first multiply by the phase factor 887 ! This allows to be left with only real operations afterwards. 888 ig=ipw1 889 890 do ipw=1,nincpw 891 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 892 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 893 do iffkg=1,nffkg 894 scalars(1,iffkg)=scalars(1,iffkg)+ar*ffkg(iffkg,ipw) 895 scalars(2,iffkg)=scalars(2,iffkg)+ai*ffkg(iffkg,ipw) 896 end do 897 ig=ig+1 898 end do 899 900 ! ******* Leaving the critical part ********************************* 901 902 ! DEBUG 903 ! write(std_out,*)' opernl2 : scalars' 904 ! do iffkg=1,10 905 ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg) 906 ! end do 907 ! stop 908 ! ENDDEBUG 909 910 if(istwf_k>=2)then 911 ! Impose parity of resulting scalar (this operation could be 912 ! replaced by direct saving of CPU time in the preceeding section) 913 do iffkg=1,nffkg 914 scalars(parity(iffkg),iffkg)=0.0d0 915 end do 916 end if 917 918 iffkg=0 919 iffkgs=nffkge+nffkgd 920 iffkgs2=nffkge+nffkgs 921 iffkgk=nffkge*2 922 do ilang=1,nlang 923 nproj=jproj(ilang) 924 if(nproj>0)then 925 ! ilang2 is the number of independent tensor components 926 ! for symmetric tensor of rank ilang-1 927 ilang2=(ilang*(ilang+1))/2 928 929 ! Loop over projectors 930 do iproj=1,nproj 931 ! Multiply by the k+G factors (tensors of various rank) 932 do ii=1,ilang2 933 ! Get the starting address for the relevant tensor 934 jj=ii+((ilang-1)*ilang*(ilang+1))/6 935 iffkg=iffkg+1 936 ! !$OMP CRITICAL (OPERNL3_1) 937 gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg) 938 gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg) 939 ! !$OMP END CRITICAL (OPERNL3_1) 940 ! Now, compute gradients, if needed. 941 if ((choice==2.or.choice==23) .and. ndgxdt==3) then 942 do mu=1,3 943 jffkg=nffkge+(iffkg-1)*3+mu 944 ! Pay attention to the use of reals and imaginary parts here ... 945 ! !$OMP CRITICAL (OPERNL3_2) 946 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 947 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 948 ! !$OMP END CRITICAL (OPERNL3_2) 949 end do 950 end if 951 if (choice==2 .and. ndgxdt==1) then 952 jffkg=nffkge+iffkg 953 ! Pay attention to the use of reals and imaginary parts here ... 954 ! !$OMP CRITICAL (OPERNL3_3) 955 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg) 956 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg) 957 ! !$OMP END CRITICAL (OPERNL3_3) 958 end if 959 if (choice==4) then 960 do mu=1,3 961 jffkg=nffkge+(iffkg-1)*9+mu 962 ! Pay attention to the use of reals and imaginary parts here ... 963 ! !$OMP CRITICAL (OPERNL3_4) 964 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 965 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 966 ! !$OMP END CRITICAL (OPERNL3_4) 967 end do 968 do mu=4,9 969 jffkg=nffkge+(iffkg-1)*9+mu 970 ! Pay attention to the use of reals and imaginary parts here ... 971 ! Also, note the multiplication by (2 pi)**2 972 ! !$OMP CRITICAL (OPERNL3_5) 973 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg) 974 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg) 975 ! !$OMP END CRITICAL (OPERNL3_5) 976 end do 977 end if 978 ! End loop on ii=1,ilang2 979 end do 980 981 if (choice==3 .or. choice==6 .or. choice==23) then 982 ! Compute additional tensors related to strain gradients 983 ! ilang3 is number of unique tensor components of rank ilang+1 984 ilang3=((ilang+2)*(ilang+3))/2 985 jjs=((ilang+1)*(ilang+2)*(ilang+3))/6 986 ! Compute strain gradient tensor components 987 do ii=1,ilang3 988 ! Note that iffkgs is also used by ddk and 2nd derivative parts 989 iffkgs=iffkgs+1 990 jj=ii+jjs 991 ! !$OMP CRITICAL (OPERNL3_6) 992 dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs) 993 dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs) 994 ! !$OMP END CRITICAL (OPERNL3_6) 995 end do 996 end if 997 998 if (choice==6) then 999 ! Compute additional tensors related to strain 2nd derivatives 1000 ! and internal strain derivatives 1001 ! ilang6 is number of unique tensor components of rank ilang+3 1002 ilang6=((ilang+4)*(ilang+5))/2 1003 jjs=((ilang+3)*(ilang+4)*(ilang+5))/6 1004 ! Compute strain gradient tensor components 1005 do ii=1,ilang6 1006 ! Note that iffkgs is also used by ddk part 1007 iffkgs2=iffkgs2+1 1008 jj=ii+jjs 1009 ! !$OMP CRITICAL (OPERNL3_6) 1010 d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2) 1011 d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2) 1012 ! !$OMP END CRITICAL (OPERNL3_6) 1013 end do 1014 1015 ! ilang4 is number of unique tensor components of rank ilang 1016 ilang4=((ilang+1)*(ilang+2))/2 1017 jjs=((ilang)*(ilang+1)*(ilang+2))/6 1018 ! Compute internal strain gradient tensor components 1019 do ii=1,ilang4 1020 iffkgs2=iffkgs2+1 1021 jj=ii+jjs 1022 ! !$OMP CRITICAL (OPERNL3_6) 1023 ! Pay attention to the use of reals and imaginary parts here ... 1024 dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2) 1025 dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2) 1026 ! !$OMP END CRITICAL (OPERNL3_6) 1027 end do 1028 1029 ! ilang5 is number of unique tensor components of rank ilang+2 1030 ilang5=((ilang+3)*(ilang+4))/2 1031 jjs=((ilang+2)*(ilang+3)*(ilang+4))/6 1032 ! Compute internal strain gradient tensor components 1033 do ii=1,ilang5 1034 iffkgs2=iffkgs2+1 1035 jj=ii+jjs 1036 ! !$OMP CRITICAL (OPERNL3_6) 1037 ! Pay attention to the use of reals and imaginary parts here ... 1038 d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2) 1039 d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2) 1040 ! !$OMP END CRITICAL (OPERNL3_6) 1041 end do 1042 end if ! choice==6 1043 1044 if (choice==5) then 1045 ! Compute additional tensors related to ddk with ffnl(:,2,...) 1046 ilangx=(ilang*(ilang+1))/2 1047 jjs=((ilang-1)*ilang*(ilang+1))/6 1048 do ii=1,ilangx 1049 ! Note that iffkgs is also used by strain part 1050 iffkgs=iffkgs+1 1051 jj=ii+jjs 1052 ! !$OMP CRITICAL (OPERNL3_7) 1053 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs) 1054 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs) 1055 ! !$OMP END CRITICAL (OPERNL3_7) 1056 end do 1057 ! Compute additional tensors related to ddk with ffnl(:,1,...) 1058 if(ilang>=2)then 1059 ilangx=((ilang-1)*ilang)/2 1060 jjs=((ilang-2)*(ilang-1)*ilang)/6 1061 do ii=1,ilangx 1062 iffkgk=iffkgk+1 1063 jj=ii+jjs 1064 ! !$OMP CRITICAL (OPERNL3_8) 1065 dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk) 1066 dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk) 1067 ! !$OMP END CRITICAL (OPERNL3_8) 1068 end do 1069 end if 1070 end if 1071 1072 ! End projector loop 1073 end do 1074 1075 ! End condition of non-zero projectors 1076 end if 1077 1078 ! End angular momentum loop 1079 end do 1080 1081 ! End loop on atoms 1082 end do 1083 1084 else if (sign==-1) then 1085 ! Application of non-local part from projected scalars 1086 ! back to reciprocal space ... 1087 ! [this section merely computes terms which add to <G|Vnl|C>; 1088 ! nothing here is needed when various gradients are being computed] 1089 1090 ! Loop on atoms 1091 do ia=1,nincat 1092 1093 ! Compute the shift eventually needed to get the phases in ph3d 1094 iaph3d=ia 1095 if(nloalg(2)>0)iaph3d=ia+ia3-1 1096 1097 ! Transfer gxa (and eventually dgxdt) in scalars with different indexing 1098 iffkg=0 1099 iffkgk=nffkge*2 1100 iffkgs=nffkge 1101 do ilang=1,nlang 1102 nproj=jproj(ilang) 1103 if (nproj>0) then 1104 ilang2=(ilang*(ilang+1))/2 1105 ilang3=((ilang+2)*(ilang+3))/2 1106 do iproj=1,nproj 1107 do ii=1,ilang2 1108 jj=ii+((ilang-1)*ilang*(ilang+1))/6 1109 iffkg=iffkg+1 1110 if(choice==1 .or. choice==3)then 1111 scalars(1,iffkg)=gxa(1,jj,ia,iproj) 1112 scalars(2,iffkg)=gxa(2,jj,ia,iproj) 1113 else if (choice==2 .and. ndgxdt==1) then 1114 jffkg=nffkge+iffkg 1115 ! Pay attention to the use of reals and imaginary parts here ... 1116 ! Also, the gxa and dgxdt arrays are switched, in order 1117 ! to give the correct combination when multiplying ffkg, 1118 ! see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]] 1119 scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj) 1120 scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj) 1121 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 1122 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 1123 else if (choice==5) then 1124 jffkg=nffkge+iffkg 1125 ! The gxa and dgxdt arrays are switched, in order 1126 ! to give the correct combination when multiplying ffkg, 1127 scalars(1,jffkg)= gxa(1,jj,ia,iproj) 1128 scalars(2,jffkg)= gxa(2,jj,ia,iproj) 1129 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 1130 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 1131 end if 1132 end do 1133 if(choice==3) then 1134 do ii=1,ilang3 1135 iffkgs=iffkgs+1 1136 jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 1137 scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj) 1138 scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj) 1139 end do 1140 end if 1141 if(ilang>=2 .and. choice==5)then 1142 do ii=1,((ilang-1)*ilang)/2 1143 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 1144 iffkgk=iffkgk+1 1145 scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj) 1146 scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj) 1147 end do 1148 end if 1149 end do 1150 end if 1151 end do 1152 1153 ! DEBUG 1154 ! if(choice==5)then 1155 ! write(std_out,*)' opernl2 : write gxa(:,...) array ' 1156 ! do ii=1,10 1157 ! write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1) 1158 ! end do 1159 ! write(std_out,*)' opernl2 : write dgxdt(:,1,...) array ' 1160 ! do ii=1,10 1161 ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1) 1162 ! end do 1163 ! write(std_out,*)' opernl2 : write dgxdt(:,2,...) array ' 1164 ! do ii=1,4 1165 ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1) 1166 ! end do 1167 ! end if 1168 1169 ! do iffkg=1,nffkg 1170 ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg) 1171 ! end do 1172 ! stop 1173 ! ENDDEBUG 1174 1175 ig=ipw1 1176 1177 ! ******* Entering the critical part ********************************* 1178 1179 do ipw=1,nincpw 1180 ar=0.0d0 1181 ai=0.0d0 1182 do iffkg=1,nffkg 1183 ar=ar+ffkg(iffkg,ipw)*scalars(1,iffkg) 1184 ai=ai+ffkg(iffkg,ipw)*scalars(2,iffkg) 1185 end do 1186 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 1187 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 1188 ig=ig+1 1189 end do 1190 1191 ! ******* Leaving the critical part ********************************* 1192 1193 ! End loop on atoms 1194 end do 1195 1196 ! End sign==-1 1197 else 1198 1199 ! Problem: sign and/or choice do not make sense 1200 write(message, '(a,2i10,a,a)' )& 1201 & ' Input sign, choice=',sign,choice,ch10,& 1202 & ' Are not compatible or allowed. ' 1203 ABI_BUG(message) 1204 end if 1205 1206 ! End loop on blocks of planewaves 1207 end do 1208 !!$OMP END DO 1209 ABI_FREE(ffkg) 1210 ABI_FREE(kpgx) 1211 ABI_FREE(parity) 1212 ABI_FREE(scalars) 1213 ABI_FREE(teffv) 1214 !!$OMP END PARALLEL 1215 1216 1217 !DEBUG 1218 !if(choice==5)then 1219 !write(std_out,*)'opernl2 : write vect(2*npw)' 1220 !do ipw=1,2 1221 !write(std_out,*)ipw,vect(1:2,ipw) 1222 !end do 1223 !write(std_out,*)'opernl2 : write ph3d' 1224 !do ipw=1,npw 1225 !write(std_out,*)ipw,ph3d(1:2,ipw,1) 1226 !end do 1227 !write(std_out,*)' opernl2 : write gxa array ' 1228 !write(std_out,*)' ang mom ,ia ' 1229 !do iproj=1,mproj 1230 !do ia=1,1 1231 !do ii=1,3 1232 !do ii=1,10 1233 !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1) 1234 !end do 1235 !end do 1236 !end do 1237 !end if 1238 !if(choice==5)then 1239 !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array ' 1240 !do ii=1,10 1241 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1) 1242 !end do 1243 !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array ' 1244 !do ii=1,4 1245 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1) 1246 !end do 1247 !stop 1248 !end if 1249 !ENDDEBUG 1250 1251 end subroutine opernl3
ABINIT/opernl4a [ Functions ]
NAME
opernl4a
FUNCTION
Operate with the non-local part of the hamiltonian, from reciprocal space to projected quantities (oprnl4b is from projected quantities to reciprocal space)
INPUTS
ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere. 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 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 mlang1 = dimensions for dgxdis1 mlang3 = one of the dimensions of the array gxa mlang4 = dimension for dgxds mlang5 = dimensions for dgxdis2 mlang6 = dimension for d2gxds2 mproj=maximum dimension for number of projection operators for each angular momentum for nonlocal pseudopotential ndgxdt=second dimension of dgxdt nffnl=third dimension of ffnl nincat = number of atoms in the subset here treated nkpg=second size of array kpg_k 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 vect(2*npw)=starting vector in reciprocal space ph3d(2,npw,matblk)=three-dimensional phase factors
OUTPUT
gxa(2,mlang3,mincat,mproj)= projected scalars if(choice==2 .or choice==4 .or. choice==5 .or. choice==23) dgxdt(2,ndgxdt,mlang3,mincat,mproj)= gradients of projected scalars wrt coords or with respect to ddk if(choice==3 .or. choice==23) dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains if(choice==6) dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars wrt coord. indexed for internal strain d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars wrt strain and coord d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars wrt strains
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
SOURCE
1323 subroutine opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 1324 & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,& 1325 & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 1326 & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,& 1327 & ntypat,ph3d,vect) 1328 1329 !Arguments ------------------------------------ 1330 !scalars 1331 integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk 1332 integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt 1333 integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat 1334 !arrays 1335 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 1336 integer,intent(in) :: nloalg(3) 1337 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 1338 real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk),vect(:,:) 1339 real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj) 1340 real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj) 1341 real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj) 1342 real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) !vz_i 1343 real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) !vz_i 1344 real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj) !vz_i 1345 1346 !Local variables------------------------------- 1347 !scalars 1348 integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2 1349 integer :: ilang3,ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs 1350 integer :: jump,mblkpw,mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2 1351 integer :: nincpw,nproj,ntens,start 1352 real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8 1353 real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8 1354 real(dp),parameter :: two_pi2=two_pi*two_pi 1355 !arrays 1356 integer,allocatable :: parity(:) 1357 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 1358 1359 ! ************************************************************************* 1360 1361 !call wrtout(std_out,"in opernl4a","COLL") 1362 1363 !mblkpw sets the size of blocks of planewaves 1364 mblkpw=NLO_MBLKPW 1365 1366 !jump governs, in fine, the use of registers in the most cpu 1367 !time consuming part of the routine. Until now, jump=8 is the maximal value. 1368 !The optimal value will be machine-dependent ! 1369 jump=4 1370 1371 !Get the actual maximum number of projectors 1372 mmproj=maxval(indlmn(3,:,itypat)) 1373 1374 !Initialisation before blocking on the plane waves 1375 1376 !Put projected scalars to zero 1377 gxa(:,:,:,1:mmproj)=0.0d0 1378 if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0 1379 if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0 1380 if (choice==6) then 1381 dgxdis(:,:,:,1:mmproj)=0.0d0 1382 d2gxdis(:,:,:,1:mmproj)=0.0d0 1383 d2gxds2(:,:,:,1:mmproj)=0.0d0 1384 end if 1385 1386 !Set up dimension of kpgx and allocate 1387 !ntens sets the maximum number of independent tensor components 1388 !over all allowed angular momenta; need 20 for spdf for tensors 1389 !up to rank 3; to handle stress tensor, need up to rank 5 1390 ntens=1 1391 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23)ntens=4 1392 if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10 1393 if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20 1394 if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35 1395 if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56 1396 if(choice==6 .and. nlang>=3)ntens=84 1397 if(choice==6 .and. nlang==4)ntens=120 1398 1399 !Set up second dimension of ffkg array, and allocate 1400 nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0 1401 do ilang=1,nlang 1402 ! Get the number of projectors for that angular momentum 1403 nproj=jproj(ilang) 1404 ! If there is a non-local part, accumulate the number of vectors needed 1405 ! The variables ilang below are the number of independent tensors of 1406 ! various ranks, the variable names being more historical than logical. 1407 ! ilang2=number of rank ilang-1 1408 ! ilang3=number of rank ilang+1 1409 ! ilang4=number of rank ilang 1410 ! ilang5=number of rank ilang+2 1411 ! ilang6=number of rank ilang+3 1412 if(nproj>0)then 1413 ilang2=(ilang*(ilang+1))/2 1414 nffkge=nffkge+nproj*ilang2 1415 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 1416 if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2 1417 if(choice==3 .or. choice==6 .or. choice==23)then 1418 ilang3=((ilang+2)*(ilang+3))/2 1419 nffkgs=nffkgs+nproj*ilang3 1420 end if 1421 if(choice==6)then 1422 ilang4=((ilang+1)*(ilang+2))/2 1423 ilang5=((ilang+3)*(ilang+4))/2 1424 ilang6=((ilang+4)*(ilang+5))/2 1425 nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6) 1426 end if 1427 end if 1428 end do 1429 nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk 1430 1431 !DEBUG 1432 !write(std_out,*)' jproj(1:nlang)',jproj(1:nlang) 1433 !write(std_out,*)' nffkg,nffkge,nffkgd,nffkgs,nffkgk',nffkg,nffkge,nffkgd,nffkgs,nffkgk 1434 !ENDDEBUG 1435 1436 !Loop on subsets of plane waves (blocking) 1437 1438 !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a 1439 !sigfault with nthreads==1 1440 !Feb 16 2012: The code does not crash anymore but it's not efficient. 1441 ! 1442 !!$OMP PARALLEL DEFAULT(PRIVATE) & 1443 !!$OMP SHARED (choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) & 1444 !!$OMP SHARED (ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) & 1445 !!$OMP SHARED (jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) & 1446 !!$OMP SHARED (mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) & 1447 !!$OMP SHARED (ntypat,ph3d,vect) & 1448 !!$OMP SHARED (mblkpw,jump,nffkgd,nffkg,nffkge,nffkgs,ntens) 1449 1450 ABI_MALLOC(ffkg,(nffkg,mblkpw)) 1451 ABI_MALLOC(parity,(nffkg)) 1452 ABI_MALLOC(kpgx,(mblkpw,ntens)) 1453 ABI_MALLOC(scalars,(2,nffkg)) 1454 ABI_MALLOC(teffv,(2,mblkpw)) 1455 1456 !!$OMP DO 1457 do ipw1=1,npw,mblkpw 1458 1459 ipw2=min(npw,ipw1+mblkpw-1) 1460 nincpw=ipw2-ipw1+1 1461 1462 ! Initialize kpgx array related to tensors defined below 1463 call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 1464 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 1465 & npw,ntens,ntypat,parity) 1466 1467 do ia=1,nincat 1468 1469 ! Compute the shift eventually needed to get the phases in ph3d 1470 iaph3d=ia 1471 if(nloalg(2)>0)iaph3d=ia+ia3-1 1472 1473 do iffkg=1,nffkg 1474 scalars(1,iffkg)=0.0d0 ; scalars(2,iffkg)=0.0d0 1475 end do 1476 1477 ! DEBUG 1478 ! write(std_out,*)'opernl4, before first time-consuming' 1479 ! write(std_out,*)'opernl4 : nffkg,nincpw=',nffkg,nincpw 1480 ! write(std_out,*)'ig,ipw,ffkg(1:4),vect(1:2)' 1481 ! ig=ipw1 1482 ! do ipw=1,nincpw 1483 ! write(std_out,'(2i4,13es11.3)' )ig,ipw,ffkg(1:min(9,nffkg),ipw),vect(1:2,ipw),ph3d(1:2,ipw,iaph3d) 1484 ! ig=ig+1 1485 ! end do 1486 ! stop 1487 ! ENDDEBUG 1488 1489 ! ******* Entering the first time-consuming part of the routine ******* 1490 1491 1492 ! First, treat small nffkg; send treat the initial phase of big 1493 ! nffkg; finally treat the loop needed for big nffkg 1494 1495 ! In the loops, first multiply by the phase factor. 1496 ! This allows to be left with only real operations afterwards. 1497 1498 ! For the time being, the maximal jump allowed is 8. 1499 1500 ! 1) Here, treat small nffkg 1501 if(nffkg<=jump)then 1502 1503 select case(nffkg) 1504 1505 case(1) 1506 1507 scr1=0.0d0 ; sci1=0.0d0 1508 ig=ipw1 1509 do ipw=1,nincpw 1510 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1511 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1512 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1513 ig=ig+1 1514 end do 1515 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1516 1517 case(2) 1518 1519 ig=ipw1 1520 scr1=0.0d0 ; sci1=0.0d0 1521 scr2=0.0d0 ; sci2=0.0d0 1522 do ipw=1,nincpw 1523 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1524 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1525 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1526 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1527 ig=ig+1 1528 end do 1529 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1530 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1531 1532 case(3) 1533 1534 ig=ipw1 1535 scr1=0.0d0 ; sci1=0.0d0 1536 scr2=0.0d0 ; sci2=0.0d0 1537 scr3=0.0d0 ; sci3=0.0d0 1538 do ipw=1,nincpw 1539 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1540 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1541 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1542 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1543 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1544 ig=ig+1 1545 end do 1546 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1547 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1548 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1549 1550 case(4) 1551 1552 ig=ipw1 1553 scr1=0.0d0 ; sci1=0.0d0 1554 scr2=0.0d0 ; sci2=0.0d0 1555 scr3=0.0d0 ; sci3=0.0d0 1556 scr4=0.0d0 ; sci4=0.0d0 1557 do ipw=1,nincpw 1558 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1559 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1560 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1561 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1562 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1563 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1564 ig=ig+1 1565 end do 1566 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1567 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1568 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1569 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1570 1571 case(5) 1572 1573 ig=ipw1 1574 scr1=0.0d0 ; sci1=0.0d0 1575 scr2=0.0d0 ; sci2=0.0d0 1576 scr3=0.0d0 ; sci3=0.0d0 1577 scr4=0.0d0 ; sci4=0.0d0 1578 scr5=0.0d0 ; sci5=0.0d0 1579 do ipw=1,nincpw 1580 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1581 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1582 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1583 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1584 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1585 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1586 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1587 ig=ig+1 1588 end do 1589 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1590 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1591 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1592 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1593 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1594 1595 case(6) 1596 1597 ig=ipw1 1598 scr1=0.0d0 ; sci1=0.0d0 1599 scr2=0.0d0 ; sci2=0.0d0 1600 scr3=0.0d0 ; sci3=0.0d0 1601 scr4=0.0d0 ; sci4=0.0d0 1602 scr5=0.0d0 ; sci5=0.0d0 1603 scr6=0.0d0 ; sci6=0.0d0 1604 do ipw=1,nincpw 1605 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1606 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1607 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1608 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1609 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1610 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1611 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1612 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 1613 ig=ig+1 1614 end do 1615 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1616 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1617 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1618 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1619 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1620 scalars(1,6)=scr6 ; scalars(2,6)=sci6 1621 1622 case(7) 1623 1624 ig=ipw1 1625 scr1=0.0d0 ; sci1=0.0d0 1626 scr2=0.0d0 ; sci2=0.0d0 1627 scr3=0.0d0 ; sci3=0.0d0 1628 scr4=0.0d0 ; sci4=0.0d0 1629 scr5=0.0d0 ; sci5=0.0d0 1630 scr6=0.0d0 ; sci6=0.0d0 1631 scr7=0.0d0 ; sci7=0.0d0 1632 do ipw=1,nincpw 1633 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1634 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1635 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1636 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1637 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1638 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1639 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1640 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 1641 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 1642 ig=ig+1 1643 end do 1644 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1645 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1646 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1647 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1648 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1649 scalars(1,6)=scr6 ; scalars(2,6)=sci6 1650 scalars(1,7)=scr7 ; scalars(2,7)=sci7 1651 1652 case(8) 1653 1654 ig=ipw1 1655 scr1=0.0d0 ; sci1=0.0d0 1656 scr2=0.0d0 ; sci2=0.0d0 1657 scr3=0.0d0 ; sci3=0.0d0 1658 scr4=0.0d0 ; sci4=0.0d0 1659 scr5=0.0d0 ; sci5=0.0d0 1660 scr6=0.0d0 ; sci6=0.0d0 1661 scr7=0.0d0 ; sci7=0.0d0 1662 scr8=0.0d0 ; sci8=0.0d0 1663 do ipw=1,nincpw 1664 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1665 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1666 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1667 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1668 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1669 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1670 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1671 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 1672 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 1673 scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw) 1674 ig=ig+1 1675 end do 1676 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1677 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1678 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1679 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1680 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1681 scalars(1,6)=scr6 ; scalars(2,6)=sci6 1682 scalars(1,7)=scr7 ; scalars(2,7)=sci7 1683 scalars(1,8)=scr8 ; scalars(2,8)=sci8 1684 1685 end select 1686 1687 else 1688 ! Now treat big nffkg 1689 1690 ! 2) Here, initialize big nffkg. The only difference with the 1691 ! preceeding case is that the intermediate results are stored. 1692 1693 select case(jump) 1694 1695 case(1) 1696 1697 scr1=0.0d0 ; sci1=0.0d0 1698 ig=ipw1 1699 do ipw=1,nincpw 1700 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1701 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1702 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1703 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1704 ig=ig+1 1705 end do 1706 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1707 1708 case(2) 1709 1710 ig=ipw1 1711 scr1=0.0d0 ; sci1=0.0d0 1712 scr2=0.0d0 ; sci2=0.0d0 1713 do ipw=1,nincpw 1714 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1715 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1716 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1717 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1718 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1719 ig=ig+1 1720 end do 1721 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1722 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1723 1724 case(3) 1725 1726 ig=ipw1 1727 scr1=0.0d0 ; sci1=0.0d0 1728 scr2=0.0d0 ; sci2=0.0d0 1729 scr3=0.0d0 ; sci3=0.0d0 1730 do ipw=1,nincpw 1731 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1732 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1733 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1734 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1735 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1736 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1737 ig=ig+1 1738 end do 1739 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1740 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1741 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1742 1743 case(4) 1744 1745 ig=ipw1 1746 scr1=0.0d0 ; sci1=0.0d0 1747 scr2=0.0d0 ; sci2=0.0d0 1748 scr3=0.0d0 ; sci3=0.0d0 1749 scr4=0.0d0 ; sci4=0.0d0 1750 do ipw=1,nincpw 1751 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1752 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1753 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1754 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1755 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1756 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1757 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1758 ig=ig+1 1759 end do 1760 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1761 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1762 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1763 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1764 1765 case(5) 1766 1767 ig=ipw1 1768 scr1=0.0d0 ; sci1=0.0d0 1769 scr2=0.0d0 ; sci2=0.0d0 1770 scr3=0.0d0 ; sci3=0.0d0 1771 scr4=0.0d0 ; sci4=0.0d0 1772 scr5=0.0d0 ; sci5=0.0d0 1773 do ipw=1,nincpw 1774 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1775 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1776 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1777 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1778 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1779 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1780 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1781 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1782 ig=ig+1 1783 end do 1784 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1785 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1786 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1787 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1788 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1789 1790 case(6) 1791 1792 ig=ipw1 1793 scr1=0.0d0 ; sci1=0.0d0 1794 scr2=0.0d0 ; sci2=0.0d0 1795 scr3=0.0d0 ; sci3=0.0d0 1796 scr4=0.0d0 ; sci4=0.0d0 1797 scr5=0.0d0 ; sci5=0.0d0 1798 scr6=0.0d0 ; sci6=0.0d0 1799 do ipw=1,nincpw 1800 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1801 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1802 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1803 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1804 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1805 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1806 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1807 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1808 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 1809 ig=ig+1 1810 end do 1811 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1812 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1813 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1814 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1815 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1816 scalars(1,6)=scr6 ; scalars(2,6)=sci6 1817 1818 case(7) 1819 1820 ig=ipw1 1821 scr1=0.0d0 ; sci1=0.0d0 1822 scr2=0.0d0 ; sci2=0.0d0 1823 scr3=0.0d0 ; sci3=0.0d0 1824 scr4=0.0d0 ; sci4=0.0d0 1825 scr5=0.0d0 ; sci5=0.0d0 1826 scr6=0.0d0 ; sci6=0.0d0 1827 scr7=0.0d0 ; sci7=0.0d0 1828 do ipw=1,nincpw 1829 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1830 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1831 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1832 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1833 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1834 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1835 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1836 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1837 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 1838 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 1839 ig=ig+1 1840 end do 1841 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1842 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1843 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1844 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1845 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1846 scalars(1,6)=scr6 ; scalars(2,6)=sci6 1847 scalars(1,7)=scr7 ; scalars(2,7)=sci7 1848 1849 case(8) 1850 1851 ig=ipw1 1852 scr1=0.0d0 ; sci1=0.0d0 1853 scr2=0.0d0 ; sci2=0.0d0 1854 scr3=0.0d0 ; sci3=0.0d0 1855 scr4=0.0d0 ; sci4=0.0d0 1856 scr5=0.0d0 ; sci5=0.0d0 1857 scr6=0.0d0 ; sci6=0.0d0 1858 scr7=0.0d0 ; sci7=0.0d0 1859 scr8=0.0d0 ; sci8=0.0d0 1860 do ipw=1,nincpw 1861 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 1862 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 1863 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 1864 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 1865 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 1866 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 1867 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 1868 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 1869 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 1870 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 1871 scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw) 1872 ig=ig+1 1873 end do 1874 scalars(1,1)=scr1 ; scalars(2,1)=sci1 1875 scalars(1,2)=scr2 ; scalars(2,2)=sci2 1876 scalars(1,3)=scr3 ; scalars(2,3)=sci3 1877 scalars(1,4)=scr4 ; scalars(2,4)=sci4 1878 scalars(1,5)=scr5 ; scalars(2,5)=sci5 1879 scalars(1,6)=scr6 ; scalars(2,6)=sci6 1880 scalars(1,7)=scr7 ; scalars(2,7)=sci7 1881 scalars(1,8)=scr8 ; scalars(2,8)=sci8 1882 1883 end select 1884 1885 ! 3) Here, do-loop for big nffkg. 1886 1887 do start=1+jump,nffkg,jump 1888 chunk=min(jump,nffkg-start+1) 1889 1890 select case(chunk) 1891 1892 case(1) 1893 1894 scr1=0.0d0 ; sci1=0.0d0 1895 do ipw=1,nincpw 1896 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 1897 scr1=scr1+ar*ffkg(start,ipw) ; sci1=sci1+ai*ffkg(start,ipw) 1898 end do 1899 scalars(1,start)=scr1 ; scalars(2,start)=sci1 1900 1901 case(2) 1902 1903 scr1=0.0d0 ; sci1=0.0d0 1904 scr2=0.0d0 ; sci2=0.0d0 1905 do ipw=1,nincpw 1906 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 1907 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 1908 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 1909 end do 1910 scalars(1,start )=scr1 ; scalars(2,start )=sci1 1911 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 1912 1913 case(3) 1914 1915 scr1=0.0d0 ; sci1=0.0d0 1916 scr2=0.0d0 ; sci2=0.0d0 1917 scr3=0.0d0 ; sci3=0.0d0 1918 do ipw=1,nincpw 1919 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 1920 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 1921 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 1922 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 1923 end do 1924 scalars(1,start )=scr1 ; scalars(2,start )=sci1 1925 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 1926 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 1927 1928 case(4) 1929 1930 scr1=0.0d0 ; sci1=0.0d0 1931 scr2=0.0d0 ; sci2=0.0d0 1932 scr3=0.0d0 ; sci3=0.0d0 1933 scr4=0.0d0 ; sci4=0.0d0 1934 do ipw=1,nincpw 1935 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 1936 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 1937 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 1938 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 1939 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 1940 end do 1941 scalars(1,start )=scr1 ; scalars(2,start )=sci1 1942 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 1943 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 1944 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 1945 1946 case(5) 1947 1948 scr1=0.0d0 ; sci1=0.0d0 1949 scr2=0.0d0 ; sci2=0.0d0 1950 scr3=0.0d0 ; sci3=0.0d0 1951 scr4=0.0d0 ; sci4=0.0d0 1952 scr5=0.0d0 ; sci5=0.0d0 1953 do ipw=1,nincpw 1954 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 1955 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 1956 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 1957 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 1958 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 1959 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 1960 end do 1961 scalars(1,start )=scr1 ; scalars(2,start )=sci1 1962 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 1963 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 1964 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 1965 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 1966 1967 case(6) 1968 1969 scr1=0.0d0 ; sci1=0.0d0 1970 scr2=0.0d0 ; sci2=0.0d0 1971 scr3=0.0d0 ; sci3=0.0d0 1972 scr4=0.0d0 ; sci4=0.0d0 1973 scr5=0.0d0 ; sci5=0.0d0 1974 scr6=0.0d0 ; sci6=0.0d0 1975 do ipw=1,nincpw 1976 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 1977 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 1978 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 1979 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 1980 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 1981 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 1982 scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw) 1983 end do 1984 scalars(1,start )=scr1 ; scalars(2,start )=sci1 1985 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 1986 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 1987 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 1988 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 1989 scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6 1990 1991 case(7) 1992 1993 scr1=0.0d0 ; sci1=0.0d0 1994 scr2=0.0d0 ; sci2=0.0d0 1995 scr3=0.0d0 ; sci3=0.0d0 1996 scr4=0.0d0 ; sci4=0.0d0 1997 scr5=0.0d0 ; sci5=0.0d0 1998 scr6=0.0d0 ; sci6=0.0d0 1999 scr7=0.0d0 ; sci7=0.0d0 2000 do ipw=1,nincpw 2001 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2002 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 2003 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 2004 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 2005 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 2006 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 2007 scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw) 2008 scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw) 2009 end do 2010 scalars(1,start )=scr1 ; scalars(2,start )=sci1 2011 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 2012 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 2013 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 2014 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 2015 scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6 2016 scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7 2017 2018 case(8) 2019 2020 scr1=0.0d0 ; sci1=0.0d0 2021 scr2=0.0d0 ; sci2=0.0d0 2022 scr3=0.0d0 ; sci3=0.0d0 2023 scr4=0.0d0 ; sci4=0.0d0 2024 scr5=0.0d0 ; sci5=0.0d0 2025 scr6=0.0d0 ; sci6=0.0d0 2026 scr7=0.0d0 ; sci7=0.0d0 2027 scr8=0.0d0 ; sci8=0.0d0 2028 do ipw=1,nincpw 2029 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2030 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 2031 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 2032 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 2033 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 2034 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 2035 scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw) 2036 scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw) 2037 scr8=scr8+ar*ffkg(start+7,ipw) ; sci8=sci8+ai*ffkg(start+7,ipw) 2038 end do 2039 scalars(1,start )=scr1 ; scalars(2,start )=sci1 2040 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 2041 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 2042 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 2043 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 2044 scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6 2045 scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7 2046 scalars(1,start+7)=scr8 ; scalars(2,start+7)=sci8 2047 2048 end select 2049 2050 ! End loop on start 2051 end do 2052 2053 ! End if statement for small or big nffkg 2054 end if 2055 2056 ! ******* Leaving the critical part ********************************* 2057 2058 ! DEBUG 2059 ! write(std_out,*)' opernl4a, write scalars ' 2060 ! do iffkg=1,nffkg 2061 ! write(std_out,*)iffkg,scalars(1:2,iffkg) 2062 ! end do 2063 ! ENDDEBUG 2064 2065 if(istwf_k>=2)then 2066 ! Impose parity of resulting scalar (this operation could be 2067 ! replaced by direct saving of CPU time in the preceeding section) 2068 do iffkg=1,nffkg 2069 scalars(parity(iffkg),iffkg)=0.0d0 2070 end do 2071 end if 2072 2073 iffkg=0 ; iffkgs=nffkge+nffkgd ; iffkgk=nffkge*2 2074 iffkgs2=nffkge+nffkgs 2075 do ilang=1,nlang 2076 nproj=jproj(ilang) 2077 if(nproj>0)then 2078 ! ilang2 is the number of independent tensor components 2079 ! for symmetric tensor of rank ilang-1 2080 ilang2=(ilang*(ilang+1))/2 2081 2082 ! Loop over projectors 2083 do iproj=1,nproj 2084 ! Multiply by the k+G factors (tensors of various rank) 2085 do ii=1,ilang2 2086 ! Get the starting address for the relevant tensor 2087 jj=ii+((ilang-1)*ilang*(ilang+1))/6 2088 iffkg=iffkg+1 2089 ! !$OMP CRITICAL (OPERNL4a_1) 2090 gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg) 2091 gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg) 2092 ! !$OMP END CRITICAL (OPERNL4a_1) 2093 ! Now, compute gradients, if needed. 2094 if ((choice==2.or.choice==23) .and. ndgxdt==3) then 2095 do mu=1,3 2096 jffkg=nffkge+(iffkg-1)*3+mu 2097 ! Pay attention to the use of reals and imaginary parts here ... 2098 ! !$OMP CRITICAL (OPERNL4a_2) 2099 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 2100 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 2101 ! !$OMP END CRITICAL (OPERNL4a_2) 2102 end do 2103 end if 2104 if (choice==2 .and. ndgxdt==1) then 2105 jffkg=nffkge+iffkg 2106 ! Pay attention to the use of reals and imaginary parts here ... 2107 ! !$OMP CRITICAL (OPERNL4a_3) 2108 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg) 2109 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg) 2110 ! !$OMP END CRITICAL (OPERNL4a_3) 2111 end if 2112 if (choice==4) then 2113 do mu=1,3 2114 jffkg=nffkge+(iffkg-1)*9+mu 2115 ! Pay attention to the use of reals and imaginary parts here ... 2116 ! !$OMP CRITICAL (OPERNL4a_4) 2117 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 2118 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 2119 ! !$OMP END CRITICAL (OPERNL4a_4) 2120 end do 2121 do mu=4,9 2122 jffkg=nffkge+(iffkg-1)*9+mu 2123 ! Pay attention to the use of reals and imaginary parts here ... 2124 ! Also, note the multiplication by (2 pi)**2 2125 ! !$OMP CRITICAL (OPERNL4a_5) 2126 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg) 2127 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg) 2128 ! !$OMP END CRITICAL (OPERNL4a_5) 2129 end do 2130 end if 2131 ! End loop on ii=1,ilang2 2132 end do 2133 2134 if ((choice==3.or.choice==23) .or. choice==6) then 2135 ! Compute additional tensors related to strain gradients 2136 ! ilang3 is number of unique tensor components of rank ilang+1 2137 ilang3=((ilang+2)*(ilang+3))/2 2138 jjs=((ilang+1)*(ilang+2)*(ilang+3))/6 2139 ! Compute strain gradient tensor components 2140 do ii=1,ilang3 2141 ! Note that iffkgs is also used by ddk and 2nd derivative parts 2142 iffkgs=iffkgs+1 2143 jj=ii+jjs 2144 ! !$OMP CRITICAL (OPERNL4a_6) 2145 dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs) 2146 dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs) 2147 ! !$OMP END CRITICAL (OPERNL4a_6) 2148 end do 2149 end if 2150 2151 if (choice==6) then 2152 ! Compute additional tensors related to strain 2nd derivatives 2153 ! and internal strain derivatives 2154 ! ilang6 is number of unique tensor components of rank ilang+3 2155 ilang6=((ilang+4)*(ilang+5))/2 2156 jjs=((ilang+3)*(ilang+4)*(ilang+5))/6 2157 ! Compute strain gradient tensor components 2158 do ii=1,ilang6 2159 ! Note that iffkgs is also used by ddk part 2160 iffkgs2=iffkgs2+1 2161 jj=ii+jjs 2162 ! !$OMP CRITICAL (OPERNL4a_7) 2163 d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2) 2164 d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2) 2165 ! !$OMP END CRITICAL (OPERNL4a_7) 2166 end do 2167 2168 ! ilang4 is number of unique tensor components of rank ilang 2169 ilang4=((ilang+1)*(ilang+2))/2 2170 jjs=((ilang)*(ilang+1)*(ilang+2))/6 2171 ! Compute internal strain gradient tensor components 2172 do ii=1,ilang4 2173 iffkgs2=iffkgs2+1 2174 jj=ii+jjs 2175 ! !$OMP CRITICAL (OPERNL4a_8) 2176 ! Pay attention to the use of reals and imaginary parts here ... 2177 dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2) 2178 dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2) 2179 ! !$OMP END CRITICAL (OPERNL4a_8) 2180 end do 2181 2182 ! ilang5 is number of unique tensor components of rank ilang+2 2183 ilang5=((ilang+3)*(ilang+4))/2 2184 jjs=((ilang+2)*(ilang+3)*(ilang+4))/6 2185 ! Compute internal strain gradient tensor components 2186 do ii=1,ilang5 2187 iffkgs2=iffkgs2+1 2188 jj=ii+jjs 2189 ! !$OMP CRITICAL (OPERNL4a_9) 2190 ! Pay attention to the use of reals and imaginary parts here ... 2191 d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2) 2192 d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2) 2193 ! !$OMP END CRITICAL (OPERNL4a_9) 2194 end do 2195 end if ! choice==6 2196 2197 if (choice==5) then 2198 ! Compute additional tensors related to ddk with ffnl(:,2,...) 2199 ilangx=(ilang*(ilang+1))/2 2200 jjs=((ilang-1)*ilang*(ilang+1))/6 2201 do ii=1,ilangx 2202 ! Note that iffkgs is also used by strain part 2203 iffkgs=iffkgs+1 2204 jj=ii+jjs 2205 ! !$OMP CRITICAL (OPERNL4a_10) 2206 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs) 2207 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs) 2208 ! !$OMP END CRITICAL (OPERNL4a_10) 2209 end do 2210 ! Compute additional tensors related to ddk with ffnl(:,1,...) 2211 if(ilang>=2)then 2212 ilangx=((ilang-1)*ilang)/2 2213 jjs=((ilang-2)*(ilang-1)*ilang)/6 2214 do ii=1,ilangx 2215 iffkgk=iffkgk+1 2216 jj=ii+jjs 2217 ! !$OMP CRITICAL (OPERNL4a_11) 2218 dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk) 2219 dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk) 2220 ! !$OMP END CRITICAL (OPERNL4a_11) 2221 end do 2222 end if 2223 end if 2224 2225 ! End projector loop 2226 end do 2227 2228 ! End condition of non-zero projectors 2229 end if 2230 2231 ! End angular momentum loop 2232 end do 2233 2234 ! End loop on atoms 2235 end do 2236 2237 ! End loop on blocks of planewaves 2238 end do 2239 !!$OMP END DO 2240 2241 ABI_FREE(ffkg) 2242 ABI_FREE(kpgx) 2243 ABI_FREE(parity) 2244 ABI_FREE(scalars) 2245 ABI_FREE(teffv) 2246 !!$OMP END PARALLEL 2247 2248 !DEBUG 2249 !write(std_out,*)' opernl4a : exit' 2250 !ENDDEBUG 2251 2252 end subroutine opernl4a
ABINIT/opernl4b [ Functions ]
NAME
opernl4b
FUNCTION
Operate with the non-local part of the hamiltonian, from projected quantities to reciprocal space
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.
SOURCE
2325 subroutine opernl4b(choice,dgxds,dgxdt,ffnl,gmet,gxa,& 2326 & ia3,idir,indlmn,ispinor,itypat,jproj,kg_k,kpg_k,kpt,& 2327 & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat,& 2328 & nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) 2329 2330 !Arguments ------------------------------------ 2331 !scalars 2332 integer,intent(in) :: choice,ia3,idir,ispinor,itypat,lmnmax,matblk !,istwf_k 2333 integer,intent(in) :: mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat !,mlang6 2334 integer,intent(in) :: nkpg,nlang,npw,ntypat 2335 !arrays 2336 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 2337 integer,intent(in) :: nloalg(3) 2338 !real(dp),intent(in) :: d2gxds2(2,mlang6,mincat,mproj) 2339 real(dp),intent(in) :: dgxds(2,mlang4,mincat,mproj) 2340 real(dp),intent(in) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) 2341 real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3) 2342 real(dp),intent(in) :: gxa(2,mlang3,mincat,mproj),kpg_k(npw,nkpg),kpt(3) 2343 real(dp),intent(in) :: ph3d(2,npw,matblk) 2344 real(dp),intent(inout) :: vect(:,:) !vz_i 2345 2346 !Local variables------------------------------- 2347 !scalars 2348 integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,ig,ii,ilang,ilang2,ilang3 2349 integer :: iproj,ipw,ipw1,ipw2,jffkg,jj,jump,mblkpw,nffkg 2350 integer :: nffkgd,nffkge,nffkgk,nffkgs,nincpw,nproj,ntens,start 2351 real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8 2352 real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8 2353 character(len=500) :: message 2354 !arrays 2355 integer,allocatable :: parity(:) 2356 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 2357 2358 ! ************************************************************************* 2359 2360 !mblkpw sets the size of blocks of planewaves 2361 mblkpw=NLO_MBLKPW 2362 2363 !jump governs, in fine, the use of registers in the most cpu 2364 !time consuming part of the routine. Until now, jump=8 is the maximal value. 2365 !The optimal value will be machine-dependent ! 2366 jump=4 2367 2368 !Initialisation before blocking on the plane waves 2369 2370 !Set up dimension of kpgx and allocate 2371 !ntens sets the maximum number of independent tensor components 2372 !over all allowed angular momenta; need 20 for spdf for tensors 2373 !up to rank 3; to handle stress tensor, need up to rank 5 2374 ntens=1 2375 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5)ntens=4 2376 if(nlang>=3 .or. choice==3)ntens=10 2377 if(nlang>=4 .or. (choice==3 .and. nlang>=2) )ntens=20 2378 if(choice==3 .and. nlang>=3)ntens=35 2379 if(choice==3 .and. nlang==4)ntens=56 2380 2381 !Set up second dimension of ffkg array, and allocate 2382 nffkg=0; nffkge=0; nffkgd=0; nffkgk=0; nffkgs=0 2383 do ilang=1,nlang 2384 ! Get the number of projectors for that angular momentum 2385 nproj=jproj(ilang) 2386 ! If there is a non-local part, accumulate the number of vectors needed 2387 if(nproj>0)then 2388 ilang2=(ilang*(ilang+1))/2 2389 nffkge=nffkge+nproj*ilang2 2390 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 2391 if(choice==2 .or. choice==4)nffkgd=nffkgd+ndgxdt*nproj*ilang2 2392 if(choice==3)then 2393 ilang3=((ilang+2)*(ilang+3))/2 2394 nffkgs=nffkgs+nproj*ilang3 2395 end if 2396 end if 2397 end do 2398 nffkg=nffkge+nffkgd+nffkgs+nffkgk 2399 2400 !!$OMP PARALLEL DEFAULT(PRIVATE) & 2401 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) & 2402 !!$OMP SHARED(itypat,jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang3,mlang4,mproj) & 2403 !!$OMP SHARED(ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) & 2404 !!$OMP SHARED(jump,nffkgd,nffkgk,nffkgs,mblkpw,nffkg,nffkge,ntens) 2405 2406 ABI_MALLOC(ffkg,(nffkg,mblkpw)) 2407 ABI_MALLOC(parity,(nffkg)) 2408 ABI_MALLOC(kpgx,(mblkpw,ntens)) 2409 ABI_MALLOC(scalars,(2,nffkg)) 2410 ABI_MALLOC(teffv,(2,mblkpw)) 2411 2412 !Loop on subsets of plane waves (blocking) 2413 !!$OMP DO 2414 do ipw1=1,npw,mblkpw 2415 2416 ipw2=min(npw,ipw1+mblkpw-1) 2417 nincpw=ipw2-ipw1+1 2418 2419 ! Initialize kpgx array related to tensors defined below 2420 call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 2421 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 2422 & npw,ntens,ntypat,parity) 2423 2424 if (choice==1 .or. choice==2 .or. choice==3 .or. choice==5) then 2425 ! Application of non-local part from projected scalars 2426 ! back to reciprocal space ... 2427 ! [this section merely computes terms which add to <G|Vnl|C>; 2428 ! nothing here is needed when various gradients are being computed] 2429 2430 ! Loop on atoms 2431 do ia=1,nincat 2432 2433 ! Compute the shift eventually needed to get the phases in ph3d 2434 iaph3d=ia 2435 if(nloalg(2)>0)iaph3d=ia+ia3-1 2436 2437 ! Transfer gxa (and eventually dgxdt) in scalars with different indexing 2438 iffkg=0 2439 iffkgk=nffkge*2 2440 iffkgs=nffkge 2441 do ilang=1,nlang 2442 nproj=jproj(ilang) 2443 if (nproj>0) then 2444 ilang2=(ilang*(ilang+1))/2 2445 ilang3=((ilang+2)*(ilang+3))/2 2446 do iproj=1,nproj 2447 do ii=1,ilang2 2448 jj=ii+((ilang-1)*ilang*(ilang+1))/6 2449 iffkg=iffkg+1 2450 if(choice==1 .or. choice==3)then 2451 scalars(1,iffkg)=gxa(1,jj,ia,iproj) 2452 scalars(2,iffkg)=gxa(2,jj,ia,iproj) 2453 else if (choice==2 .and. ndgxdt==1) then 2454 jffkg=nffkge+iffkg 2455 ! Pay attention to the use of reals and imaginary parts here ... 2456 ! Also, the gxa and dgxdt arrays are switched, in order 2457 ! to give the correct combination when multiplying ffkg, 2458 ! see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]] 2459 scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj) 2460 scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj) 2461 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 2462 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 2463 else if (choice==5) then 2464 jffkg=nffkge+iffkg 2465 ! The gxa and dgxdt arrays are switched, in order 2466 ! to give the correct combination when multiplying ffkg, 2467 scalars(1,jffkg)= gxa(1,jj,ia,iproj) 2468 scalars(2,jffkg)= gxa(2,jj,ia,iproj) 2469 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 2470 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 2471 end if 2472 end do 2473 if(choice==3) then 2474 do ii=1,ilang3 2475 iffkgs=iffkgs+1 2476 jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 2477 scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj) 2478 scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj) 2479 end do 2480 end if 2481 if(ilang>=2 .and. choice==5)then 2482 do ii=1,((ilang-1)*ilang)/2 2483 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 2484 iffkgk=iffkgk+1 2485 scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj) 2486 scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj) 2487 end do 2488 end if 2489 end do 2490 end if 2491 end do 2492 2493 ! DEBUG 2494 ! write(std_out,*)' opernl4b, write scalars ' 2495 ! do iffkg=1,nffkg 2496 ! write(std_out,*)iffkg,scalars(1:2,iffkg) 2497 ! end do 2498 ! ENDDEBUG 2499 2500 ! ******* Entering the second critical part **************************** 2501 2502 ! First, treat small nffkg; send treat the loop needed for big nffkg; 2503 ! finally treat the end of the loop needed for big nffkg 2504 2505 ! For the time being, the maximal jump allowed is 8. 2506 2507 ! 1) Here, treat small nffkg 2508 if(nffkg<=jump)then 2509 2510 select case(nffkg) 2511 2512 case(1) 2513 2514 scr1=scalars(1,1) ; sci1=scalars(2,1) 2515 ig=ipw1 2516 do ipw=1,nincpw 2517 ar=ffkg(1,ipw)*scr1 ; ai=ffkg(1,ipw)*sci1 2518 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2519 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2520 ig=ig+1 2521 end do 2522 2523 case(2) 2524 2525 scr1=scalars(1,1) ; sci1=scalars(2,1) 2526 scr2=scalars(1,2) ; sci2=scalars(2,2) 2527 ig=ipw1 2528 do ipw=1,nincpw 2529 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2530 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2531 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2532 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2533 ig=ig+1 2534 end do 2535 2536 case(3) 2537 2538 scr1=scalars(1,1) ; sci1=scalars(2,1) 2539 scr2=scalars(1,2) ; sci2=scalars(2,2) 2540 scr3=scalars(1,3) ; sci3=scalars(2,3) 2541 ig=ipw1 2542 do ipw=1,nincpw 2543 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2544 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2545 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 2546 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2547 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2548 ig=ig+1 2549 end do 2550 2551 case(4) 2552 2553 scr1=scalars(1,1) ; sci1=scalars(2,1) 2554 scr2=scalars(1,2) ; sci2=scalars(2,2) 2555 scr3=scalars(1,3) ; sci3=scalars(2,3) 2556 scr4=scalars(1,4) ; sci4=scalars(2,4) 2557 ig=ipw1 2558 do ipw=1,nincpw 2559 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2560 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2561 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 2562 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 2563 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2564 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2565 ig=ig+1 2566 end do 2567 2568 case(5) 2569 2570 scr1=scalars(1,1) ; sci1=scalars(2,1) 2571 scr2=scalars(1,2) ; sci2=scalars(2,2) 2572 scr3=scalars(1,3) ; sci3=scalars(2,3) 2573 scr4=scalars(1,4) ; sci4=scalars(2,4) 2574 scr5=scalars(1,5) ; sci5=scalars(2,5) 2575 ig=ipw1 2576 do ipw=1,nincpw 2577 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2578 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2579 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 2580 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 2581 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 2582 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2583 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2584 ig=ig+1 2585 end do 2586 2587 case(6) 2588 2589 scr1=scalars(1,1) ; sci1=scalars(2,1) 2590 scr2=scalars(1,2) ; sci2=scalars(2,2) 2591 scr3=scalars(1,3) ; sci3=scalars(2,3) 2592 scr4=scalars(1,4) ; sci4=scalars(2,4) 2593 scr5=scalars(1,5) ; sci5=scalars(2,5) 2594 scr6=scalars(1,6) ; sci6=scalars(2,6) 2595 ig=ipw1 2596 do ipw=1,nincpw 2597 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2598 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2599 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 2600 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 2601 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 2602 ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6 2603 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2604 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2605 ig=ig+1 2606 end do 2607 2608 case(7) 2609 2610 scr1=scalars(1,1) ; sci1=scalars(2,1) 2611 scr2=scalars(1,2) ; sci2=scalars(2,2) 2612 scr3=scalars(1,3) ; sci3=scalars(2,3) 2613 scr4=scalars(1,4) ; sci4=scalars(2,4) 2614 scr5=scalars(1,5) ; sci5=scalars(2,5) 2615 scr6=scalars(1,6) ; sci6=scalars(2,6) 2616 scr7=scalars(1,7) ; sci7=scalars(2,7) 2617 ig=ipw1 2618 do ipw=1,nincpw 2619 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2620 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2621 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 2622 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 2623 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 2624 ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6 2625 ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7 2626 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2627 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2628 ig=ig+1 2629 end do 2630 2631 case(8) 2632 2633 scr1=scalars(1,1) ; sci1=scalars(2,1) 2634 scr2=scalars(1,2) ; sci2=scalars(2,2) 2635 scr3=scalars(1,3) ; sci3=scalars(2,3) 2636 scr4=scalars(1,4) ; sci4=scalars(2,4) 2637 scr5=scalars(1,5) ; sci5=scalars(2,5) 2638 scr6=scalars(1,6) ; sci6=scalars(2,6) 2639 scr7=scalars(1,7) ; sci7=scalars(2,7) 2640 scr8=scalars(1,8) ; sci8=scalars(2,8) 2641 ig=ipw1 2642 do ipw=1,nincpw 2643 ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1 2644 ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2 2645 ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3 2646 ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4 2647 ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5 2648 ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6 2649 ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7 2650 ar=ar+ffkg(8,ipw)*scr8 ; ai=ai+ffkg(8,ipw)*sci8 2651 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2652 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2653 ig=ig+1 2654 end do 2655 2656 2657 end select 2658 2659 else 2660 2661 do ipw=1,nincpw 2662 teffv(1,ipw)=0.0d0 ; teffv(2,ipw)=0.0d0 2663 end do 2664 2665 ! 2) Here treart the loop for big nffkg 2666 do start=1,nffkg-jump,jump 2667 chunk=min(jump,nffkg-jump-start+1) 2668 2669 select case(chunk) 2670 2671 case(1) 2672 2673 scr1=scalars(1,start) ; sci1=scalars(2,start) 2674 do ipw=1,nincpw 2675 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2676 ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1 2677 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2678 end do 2679 2680 case(2) 2681 2682 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2683 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2684 do ipw=1,nincpw 2685 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2686 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2687 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2688 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2689 end do 2690 2691 case(3) 2692 2693 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2694 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2695 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2696 do ipw=1,nincpw 2697 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2698 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2699 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2700 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2701 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2702 end do 2703 2704 case(4) 2705 2706 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2707 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2708 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2709 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2710 do ipw=1,nincpw 2711 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2712 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2713 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2714 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2715 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2716 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2717 end do 2718 2719 case(5) 2720 2721 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2722 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2723 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2724 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2725 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2726 do ipw=1,nincpw 2727 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2728 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2729 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2730 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2731 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2732 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2733 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2734 end do 2735 2736 case(6) 2737 2738 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2739 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2740 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2741 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2742 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2743 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 2744 do ipw=1,nincpw 2745 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2746 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2747 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2748 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2749 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2750 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2751 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 2752 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2753 end do 2754 2755 case(7) 2756 2757 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2758 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2759 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2760 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2761 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2762 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 2763 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 2764 do ipw=1,nincpw 2765 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2766 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2767 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2768 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2769 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2770 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2771 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 2772 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 2773 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2774 end do 2775 2776 case(8) 2777 2778 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2779 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2780 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2781 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2782 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2783 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 2784 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 2785 scr8=scalars(1,start+7) ; sci8=scalars(2,start+7) 2786 do ipw=1,nincpw 2787 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2788 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2789 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2790 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2791 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2792 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2793 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 2794 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 2795 ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8 2796 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 2797 end do 2798 2799 end select 2800 2801 end do 2802 2803 ! 3) Treat the end of the loops 2804 2805 start=nffkg-jump+1 2806 2807 select case(jump) 2808 2809 case(1) 2810 2811 scr1=scalars(1,start) ; sci1=scalars(2,start) 2812 ig=ipw1 2813 do ipw=1,nincpw 2814 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2815 ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1 2816 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2817 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2818 ig=ig+1 2819 end do 2820 2821 case(2) 2822 2823 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2824 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2825 ig=ipw1 2826 do ipw=1,nincpw 2827 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2828 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2829 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2830 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2831 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2832 ig=ig+1 2833 end do 2834 2835 case(3) 2836 2837 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2838 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2839 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2840 ig=ipw1 2841 do ipw=1,nincpw 2842 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2843 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2844 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2845 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2846 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2847 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2848 ig=ig+1 2849 end do 2850 2851 case(4) 2852 2853 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2854 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2855 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2856 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2857 ig=ipw1 2858 do ipw=1,nincpw 2859 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2860 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2861 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2862 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2863 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2864 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2865 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2866 ig=ig+1 2867 end do 2868 2869 case(5) 2870 2871 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2872 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2873 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2874 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2875 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2876 ig=ipw1 2877 do ipw=1,nincpw 2878 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2879 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2880 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2881 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2882 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2883 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2884 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2885 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2886 ig=ig+1 2887 end do 2888 2889 case(6) 2890 2891 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2892 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2893 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2894 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2895 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2896 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 2897 ig=ipw1 2898 do ipw=1,nincpw 2899 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2900 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2901 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2902 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2903 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2904 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2905 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 2906 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2907 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2908 ig=ig+1 2909 end do 2910 2911 case(7) 2912 2913 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2914 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2915 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2916 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2917 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2918 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 2919 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 2920 ig=ipw1 2921 do ipw=1,nincpw 2922 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2923 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2924 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2925 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2926 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2927 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2928 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 2929 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 2930 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2931 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2932 ig=ig+1 2933 end do 2934 2935 case(8) 2936 2937 scr1=scalars(1,start ) ; sci1=scalars(2,start ) 2938 scr2=scalars(1,start+1) ; sci2=scalars(2,start+1) 2939 scr3=scalars(1,start+2) ; sci3=scalars(2,start+2) 2940 scr4=scalars(1,start+3) ; sci4=scalars(2,start+3) 2941 scr5=scalars(1,start+4) ; sci5=scalars(2,start+4) 2942 scr6=scalars(1,start+5) ; sci6=scalars(2,start+5) 2943 scr7=scalars(1,start+6) ; sci7=scalars(2,start+6) 2944 scr8=scalars(1,start+7) ; sci8=scalars(2,start+7) 2945 ig=ipw1 2946 do ipw=1,nincpw 2947 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 2948 ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1 2949 ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2 2950 ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3 2951 ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4 2952 ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5 2953 ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6 2954 ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7 2955 ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8 2956 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 2957 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 2958 ig=ig+1 2959 end do 2960 2961 end select 2962 2963 ! End if statement for small or big nffkg 2964 end if 2965 2966 ! ******* Leaving the critical part ********************************* 2967 2968 ! End loop on atoms 2969 end do 2970 2971 ! End choice==1 or choice==2 or choice==3 2972 else 2973 ! Problem: choice does not make sense 2974 write(message,'(a,i0,a)' )' Input choice=',choice,' not allowed. ' 2975 ABI_BUG(message) 2976 end if 2977 2978 ! End loop on blocks of planewaves 2979 end do 2980 !!$OMP END DO 2981 2982 ABI_FREE(ffkg) 2983 ABI_FREE(kpgx) 2984 ABI_FREE(parity) 2985 ABI_FREE(scalars) 2986 ABI_FREE(teffv) 2987 !!$OMP END PARALLEL 2988 2989 end subroutine opernl4b