TABLE OF CONTENTS
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)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
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
PARENTS
nonlop_pl
CHILDREN
dfpt_mkffkg
SOURCE
83 #if defined HAVE_CONFIG_H 84 #include "config.h" 85 #endif 86 87 #include "abi_common.h" 88 89 subroutine opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 90 & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,& 91 & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 92 & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,& 93 & ntypat,ph3d,vect) 94 95 use defs_basis 96 use m_profiling_abi 97 98 !This section has been created automatically by the script Abilint (TD). 99 !Do not modify the following lines by hand. 100 #undef ABI_FUNC 101 #define ABI_FUNC 'opernl4a' 102 use interfaces_66_nonlocal, except_this_one => opernl4a 103 !End of the abilint section 104 105 implicit none 106 107 !Arguments ------------------------------------ 108 !scalars 109 integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk 110 integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt 111 integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat 112 !arrays 113 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 114 integer,intent(in) :: nloalg(3) 115 real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg) 116 real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk),vect(:,:) 117 real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj) 118 real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj) 119 real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj) 120 real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) !vz_i 121 real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) !vz_i 122 real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj) !vz_i 123 124 !Local variables------------------------------- 125 !scalars 126 integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2 127 integer :: ilang3,ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs 128 integer :: jump,mblkpw,mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2 129 integer :: nincpw,nproj,ntens,start 130 real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8 131 real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8 132 real(dp),parameter :: two_pi2=two_pi*two_pi 133 !arrays 134 integer,allocatable :: parity(:) 135 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 136 137 ! ************************************************************************* 138 139 !call wrtout(std_out,"in opernl4a","COLL") 140 141 !mblkpw sets the size of blocks of planewaves 142 mblkpw=NLO_MBLKPW 143 144 !jump governs, in fine, the use of registers in the most cpu 145 !time consuming part of the routine. Until now, jump=8 is the maximal value. 146 !The optimal value will be machine-dependent ! 147 jump=4 148 149 !Get the actual maximum number of projectors 150 mmproj=maxval(indlmn(3,:,itypat)) 151 152 !Initialisation before blocking on the plane waves 153 154 !Put projected scalars to zero 155 gxa(:,:,:,1:mmproj)=0.0d0 156 if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0 157 if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0 158 if (choice==6) then 159 dgxdis(:,:,:,1:mmproj)=0.0d0 160 d2gxdis(:,:,:,1:mmproj)=0.0d0 161 d2gxds2(:,:,:,1:mmproj)=0.0d0 162 end if 163 164 !Set up dimension of kpgx and allocate 165 !ntens sets the maximum number of independent tensor components 166 !over all allowed angular momenta; need 20 for spdf for tensors 167 !up to rank 3; to handle stress tensor, need up to rank 5 168 ntens=1 169 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23)ntens=4 170 if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10 171 if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20 172 if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35 173 if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56 174 if(choice==6 .and. nlang>=3)ntens=84 175 if(choice==6 .and. nlang==4)ntens=120 176 177 !Set up second dimension of ffkg array, and allocate 178 nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0 179 do ilang=1,nlang 180 ! Get the number of projectors for that angular momentum 181 nproj=jproj(ilang) 182 ! If there is a non-local part, accumulate the number of vectors needed 183 ! The variables ilang below are the number of independent tensors of 184 ! various ranks, the variable names being more historical than logical. 185 ! ilang2=number of rank ilang-1 186 ! ilang3=number of rank ilang+1 187 ! ilang4=number of rank ilang 188 ! ilang5=number of rank ilang+2 189 ! ilang6=number of rank ilang+3 190 if(nproj>0)then 191 ilang2=(ilang*(ilang+1))/2 192 nffkge=nffkge+nproj*ilang2 193 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 194 if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2 195 if(choice==3 .or. choice==6 .or. choice==23)then 196 ilang3=((ilang+2)*(ilang+3))/2 197 nffkgs=nffkgs+nproj*ilang3 198 end if 199 if(choice==6)then 200 ilang4=((ilang+1)*(ilang+2))/2 201 ilang5=((ilang+3)*(ilang+4))/2 202 ilang6=((ilang+4)*(ilang+5))/2 203 nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6) 204 end if 205 end if 206 end do 207 nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk 208 209 !DEBUG 210 !write(std_out,*)' jproj(1:nlang)',jproj(1:nlang) 211 !write(std_out,*)' nffkg,nffkge,nffkgd,nffkgs,nffkgk',nffkg,nffkge,nffkgd,nffkgs,nffkgk 212 !ENDDEBUG 213 214 !Loop on subsets of plane waves (blocking) 215 216 !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a 217 !sigfault with nthreads==1 218 !Feb 16 2012: The code does not crash anymore but it's not efficient. 219 ! 220 !!$OMP PARALLEL DEFAULT(PRIVATE) & 221 !!$OMP SHARED (choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) & 222 !!$OMP SHARED (ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) & 223 !!$OMP SHARED (jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) & 224 !!$OMP SHARED (mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) & 225 !!$OMP SHARED (ntypat,ph3d,vect) & 226 !!$OMP SHARED (mblkpw,jump,nffkgd,nffkg,nffkge,nffkgs,ntens) 227 228 ABI_ALLOCATE(ffkg,(nffkg,mblkpw)) 229 ABI_ALLOCATE(parity,(nffkg)) 230 ABI_ALLOCATE(kpgx,(mblkpw,ntens)) 231 ABI_ALLOCATE(scalars,(2,nffkg)) 232 ABI_ALLOCATE(teffv,(2,mblkpw)) 233 234 !!$OMP DO 235 do ipw1=1,npw,mblkpw 236 237 ipw2=min(npw,ipw1+mblkpw-1) 238 nincpw=ipw2-ipw1+1 239 240 ! Initialize kpgx array related to tensors defined below 241 call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 242 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 243 & npw,ntens,ntypat,parity) 244 245 do ia=1,nincat 246 247 ! Compute the shift eventually needed to get the phases in ph3d 248 iaph3d=ia 249 if(nloalg(2)>0)iaph3d=ia+ia3-1 250 251 do iffkg=1,nffkg 252 scalars(1,iffkg)=0.0d0 ; scalars(2,iffkg)=0.0d0 253 end do 254 255 ! DEBUG 256 ! write(std_out,*)'opernl4, before first time-consuming' 257 ! write(std_out,*)'opernl4 : nffkg,nincpw=',nffkg,nincpw 258 ! write(std_out,*)'ig,ipw,ffkg(1:4),vect(1:2)' 259 ! ig=ipw1 260 ! do ipw=1,nincpw 261 ! write(std_out,'(2i4,13es11.3)' )ig,ipw,ffkg(1:min(9,nffkg),ipw),vect(1:2,ipw),ph3d(1:2,ipw,iaph3d) 262 ! ig=ig+1 263 ! end do 264 ! stop 265 ! ENDDEBUG 266 267 ! ******* Entering the first time-consuming part of the routine ******* 268 269 270 ! First, treat small nffkg; send treat the initial phase of big 271 ! nffkg; finally treat the loop needed for big nffkg 272 273 ! In the loops, first multiply by the phase factor. 274 ! This allows to be left with only real operations afterwards. 275 276 ! For the time being, the maximal jump allowed is 8. 277 278 ! 1) Here, treat small nffkg 279 if(nffkg<=jump)then 280 281 select case(nffkg) 282 283 case(1) 284 285 scr1=0.0d0 ; sci1=0.0d0 286 ig=ipw1 287 do ipw=1,nincpw 288 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 289 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 290 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 291 ig=ig+1 292 end do 293 scalars(1,1)=scr1 ; scalars(2,1)=sci1 294 295 case(2) 296 297 ig=ipw1 298 scr1=0.0d0 ; sci1=0.0d0 299 scr2=0.0d0 ; sci2=0.0d0 300 do ipw=1,nincpw 301 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 302 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 303 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 304 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 305 ig=ig+1 306 end do 307 scalars(1,1)=scr1 ; scalars(2,1)=sci1 308 scalars(1,2)=scr2 ; scalars(2,2)=sci2 309 310 case(3) 311 312 ig=ipw1 313 scr1=0.0d0 ; sci1=0.0d0 314 scr2=0.0d0 ; sci2=0.0d0 315 scr3=0.0d0 ; sci3=0.0d0 316 do ipw=1,nincpw 317 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 318 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 319 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 320 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 321 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 322 ig=ig+1 323 end do 324 scalars(1,1)=scr1 ; scalars(2,1)=sci1 325 scalars(1,2)=scr2 ; scalars(2,2)=sci2 326 scalars(1,3)=scr3 ; scalars(2,3)=sci3 327 328 case(4) 329 330 ig=ipw1 331 scr1=0.0d0 ; sci1=0.0d0 332 scr2=0.0d0 ; sci2=0.0d0 333 scr3=0.0d0 ; sci3=0.0d0 334 scr4=0.0d0 ; sci4=0.0d0 335 do ipw=1,nincpw 336 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 337 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 338 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 339 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 340 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 341 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 342 ig=ig+1 343 end do 344 scalars(1,1)=scr1 ; scalars(2,1)=sci1 345 scalars(1,2)=scr2 ; scalars(2,2)=sci2 346 scalars(1,3)=scr3 ; scalars(2,3)=sci3 347 scalars(1,4)=scr4 ; scalars(2,4)=sci4 348 349 case(5) 350 351 ig=ipw1 352 scr1=0.0d0 ; sci1=0.0d0 353 scr2=0.0d0 ; sci2=0.0d0 354 scr3=0.0d0 ; sci3=0.0d0 355 scr4=0.0d0 ; sci4=0.0d0 356 scr5=0.0d0 ; sci5=0.0d0 357 do ipw=1,nincpw 358 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 359 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 360 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 361 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 362 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 363 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 364 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 365 ig=ig+1 366 end do 367 scalars(1,1)=scr1 ; scalars(2,1)=sci1 368 scalars(1,2)=scr2 ; scalars(2,2)=sci2 369 scalars(1,3)=scr3 ; scalars(2,3)=sci3 370 scalars(1,4)=scr4 ; scalars(2,4)=sci4 371 scalars(1,5)=scr5 ; scalars(2,5)=sci5 372 373 case(6) 374 375 ig=ipw1 376 scr1=0.0d0 ; sci1=0.0d0 377 scr2=0.0d0 ; sci2=0.0d0 378 scr3=0.0d0 ; sci3=0.0d0 379 scr4=0.0d0 ; sci4=0.0d0 380 scr5=0.0d0 ; sci5=0.0d0 381 scr6=0.0d0 ; sci6=0.0d0 382 do ipw=1,nincpw 383 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 384 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 385 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 386 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 387 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 388 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 389 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 390 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 391 ig=ig+1 392 end do 393 scalars(1,1)=scr1 ; scalars(2,1)=sci1 394 scalars(1,2)=scr2 ; scalars(2,2)=sci2 395 scalars(1,3)=scr3 ; scalars(2,3)=sci3 396 scalars(1,4)=scr4 ; scalars(2,4)=sci4 397 scalars(1,5)=scr5 ; scalars(2,5)=sci5 398 scalars(1,6)=scr6 ; scalars(2,6)=sci6 399 400 case(7) 401 402 ig=ipw1 403 scr1=0.0d0 ; sci1=0.0d0 404 scr2=0.0d0 ; sci2=0.0d0 405 scr3=0.0d0 ; sci3=0.0d0 406 scr4=0.0d0 ; sci4=0.0d0 407 scr5=0.0d0 ; sci5=0.0d0 408 scr6=0.0d0 ; sci6=0.0d0 409 scr7=0.0d0 ; sci7=0.0d0 410 do ipw=1,nincpw 411 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 412 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 413 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 414 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 415 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 416 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 417 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 418 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 419 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 420 ig=ig+1 421 end do 422 scalars(1,1)=scr1 ; scalars(2,1)=sci1 423 scalars(1,2)=scr2 ; scalars(2,2)=sci2 424 scalars(1,3)=scr3 ; scalars(2,3)=sci3 425 scalars(1,4)=scr4 ; scalars(2,4)=sci4 426 scalars(1,5)=scr5 ; scalars(2,5)=sci5 427 scalars(1,6)=scr6 ; scalars(2,6)=sci6 428 scalars(1,7)=scr7 ; scalars(2,7)=sci7 429 430 case(8) 431 432 ig=ipw1 433 scr1=0.0d0 ; sci1=0.0d0 434 scr2=0.0d0 ; sci2=0.0d0 435 scr3=0.0d0 ; sci3=0.0d0 436 scr4=0.0d0 ; sci4=0.0d0 437 scr5=0.0d0 ; sci5=0.0d0 438 scr6=0.0d0 ; sci6=0.0d0 439 scr7=0.0d0 ; sci7=0.0d0 440 scr8=0.0d0 ; sci8=0.0d0 441 do ipw=1,nincpw 442 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 443 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 444 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 445 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 446 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 447 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 448 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 449 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 450 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 451 scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw) 452 ig=ig+1 453 end do 454 scalars(1,1)=scr1 ; scalars(2,1)=sci1 455 scalars(1,2)=scr2 ; scalars(2,2)=sci2 456 scalars(1,3)=scr3 ; scalars(2,3)=sci3 457 scalars(1,4)=scr4 ; scalars(2,4)=sci4 458 scalars(1,5)=scr5 ; scalars(2,5)=sci5 459 scalars(1,6)=scr6 ; scalars(2,6)=sci6 460 scalars(1,7)=scr7 ; scalars(2,7)=sci7 461 scalars(1,8)=scr8 ; scalars(2,8)=sci8 462 463 end select 464 465 else 466 ! Now treat big nffkg 467 468 ! 2) Here, initialize big nffkg. The only difference with the 469 ! preceeding case is that the intermediate results are stored. 470 471 select case(jump) 472 473 case(1) 474 475 scr1=0.0d0 ; sci1=0.0d0 476 ig=ipw1 477 do ipw=1,nincpw 478 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 479 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 480 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 481 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 482 ig=ig+1 483 end do 484 scalars(1,1)=scr1 ; scalars(2,1)=sci1 485 486 case(2) 487 488 ig=ipw1 489 scr1=0.0d0 ; sci1=0.0d0 490 scr2=0.0d0 ; sci2=0.0d0 491 do ipw=1,nincpw 492 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 493 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 494 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 495 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 496 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 497 ig=ig+1 498 end do 499 scalars(1,1)=scr1 ; scalars(2,1)=sci1 500 scalars(1,2)=scr2 ; scalars(2,2)=sci2 501 502 case(3) 503 504 ig=ipw1 505 scr1=0.0d0 ; sci1=0.0d0 506 scr2=0.0d0 ; sci2=0.0d0 507 scr3=0.0d0 ; sci3=0.0d0 508 do ipw=1,nincpw 509 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 510 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 511 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 512 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 513 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 514 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 515 ig=ig+1 516 end do 517 scalars(1,1)=scr1 ; scalars(2,1)=sci1 518 scalars(1,2)=scr2 ; scalars(2,2)=sci2 519 scalars(1,3)=scr3 ; scalars(2,3)=sci3 520 521 case(4) 522 523 ig=ipw1 524 scr1=0.0d0 ; sci1=0.0d0 525 scr2=0.0d0 ; sci2=0.0d0 526 scr3=0.0d0 ; sci3=0.0d0 527 scr4=0.0d0 ; sci4=0.0d0 528 do ipw=1,nincpw 529 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 530 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 531 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 532 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 533 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 534 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 535 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 536 ig=ig+1 537 end do 538 scalars(1,1)=scr1 ; scalars(2,1)=sci1 539 scalars(1,2)=scr2 ; scalars(2,2)=sci2 540 scalars(1,3)=scr3 ; scalars(2,3)=sci3 541 scalars(1,4)=scr4 ; scalars(2,4)=sci4 542 543 case(5) 544 545 ig=ipw1 546 scr1=0.0d0 ; sci1=0.0d0 547 scr2=0.0d0 ; sci2=0.0d0 548 scr3=0.0d0 ; sci3=0.0d0 549 scr4=0.0d0 ; sci4=0.0d0 550 scr5=0.0d0 ; sci5=0.0d0 551 do ipw=1,nincpw 552 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 553 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 554 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 555 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 556 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 557 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 558 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 559 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 560 ig=ig+1 561 end do 562 scalars(1,1)=scr1 ; scalars(2,1)=sci1 563 scalars(1,2)=scr2 ; scalars(2,2)=sci2 564 scalars(1,3)=scr3 ; scalars(2,3)=sci3 565 scalars(1,4)=scr4 ; scalars(2,4)=sci4 566 scalars(1,5)=scr5 ; scalars(2,5)=sci5 567 568 case(6) 569 570 ig=ipw1 571 scr1=0.0d0 ; sci1=0.0d0 572 scr2=0.0d0 ; sci2=0.0d0 573 scr3=0.0d0 ; sci3=0.0d0 574 scr4=0.0d0 ; sci4=0.0d0 575 scr5=0.0d0 ; sci5=0.0d0 576 scr6=0.0d0 ; sci6=0.0d0 577 do ipw=1,nincpw 578 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 579 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 580 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 581 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 582 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 583 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 584 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 585 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 586 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 587 ig=ig+1 588 end do 589 scalars(1,1)=scr1 ; scalars(2,1)=sci1 590 scalars(1,2)=scr2 ; scalars(2,2)=sci2 591 scalars(1,3)=scr3 ; scalars(2,3)=sci3 592 scalars(1,4)=scr4 ; scalars(2,4)=sci4 593 scalars(1,5)=scr5 ; scalars(2,5)=sci5 594 scalars(1,6)=scr6 ; scalars(2,6)=sci6 595 596 case(7) 597 598 ig=ipw1 599 scr1=0.0d0 ; sci1=0.0d0 600 scr2=0.0d0 ; sci2=0.0d0 601 scr3=0.0d0 ; sci3=0.0d0 602 scr4=0.0d0 ; sci4=0.0d0 603 scr5=0.0d0 ; sci5=0.0d0 604 scr6=0.0d0 ; sci6=0.0d0 605 scr7=0.0d0 ; sci7=0.0d0 606 do ipw=1,nincpw 607 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 608 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 609 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 610 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 611 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 612 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 613 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 614 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 615 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 616 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 617 ig=ig+1 618 end do 619 scalars(1,1)=scr1 ; scalars(2,1)=sci1 620 scalars(1,2)=scr2 ; scalars(2,2)=sci2 621 scalars(1,3)=scr3 ; scalars(2,3)=sci3 622 scalars(1,4)=scr4 ; scalars(2,4)=sci4 623 scalars(1,5)=scr5 ; scalars(2,5)=sci5 624 scalars(1,6)=scr6 ; scalars(2,6)=sci6 625 scalars(1,7)=scr7 ; scalars(2,7)=sci7 626 627 case(8) 628 629 ig=ipw1 630 scr1=0.0d0 ; sci1=0.0d0 631 scr2=0.0d0 ; sci2=0.0d0 632 scr3=0.0d0 ; sci3=0.0d0 633 scr4=0.0d0 ; sci4=0.0d0 634 scr5=0.0d0 ; sci5=0.0d0 635 scr6=0.0d0 ; sci6=0.0d0 636 scr7=0.0d0 ; sci7=0.0d0 637 scr8=0.0d0 ; sci8=0.0d0 638 do ipw=1,nincpw 639 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 640 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 641 teffv(1,ipw)=ar ; teffv(2,ipw)=ai 642 scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw) 643 scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw) 644 scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw) 645 scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw) 646 scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw) 647 scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw) 648 scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw) 649 scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw) 650 ig=ig+1 651 end do 652 scalars(1,1)=scr1 ; scalars(2,1)=sci1 653 scalars(1,2)=scr2 ; scalars(2,2)=sci2 654 scalars(1,3)=scr3 ; scalars(2,3)=sci3 655 scalars(1,4)=scr4 ; scalars(2,4)=sci4 656 scalars(1,5)=scr5 ; scalars(2,5)=sci5 657 scalars(1,6)=scr6 ; scalars(2,6)=sci6 658 scalars(1,7)=scr7 ; scalars(2,7)=sci7 659 scalars(1,8)=scr8 ; scalars(2,8)=sci8 660 661 end select 662 663 ! 3) Here, do-loop for big nffkg. 664 665 do start=1+jump,nffkg,jump 666 chunk=min(jump,nffkg-start+1) 667 668 select case(chunk) 669 670 case(1) 671 672 scr1=0.0d0 ; sci1=0.0d0 673 do ipw=1,nincpw 674 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 675 scr1=scr1+ar*ffkg(start,ipw) ; sci1=sci1+ai*ffkg(start,ipw) 676 end do 677 scalars(1,start)=scr1 ; scalars(2,start)=sci1 678 679 case(2) 680 681 scr1=0.0d0 ; sci1=0.0d0 682 scr2=0.0d0 ; sci2=0.0d0 683 do ipw=1,nincpw 684 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 685 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 686 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 687 end do 688 scalars(1,start )=scr1 ; scalars(2,start )=sci1 689 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 690 691 case(3) 692 693 scr1=0.0d0 ; sci1=0.0d0 694 scr2=0.0d0 ; sci2=0.0d0 695 scr3=0.0d0 ; sci3=0.0d0 696 do ipw=1,nincpw 697 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 698 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 699 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 700 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 701 end do 702 scalars(1,start )=scr1 ; scalars(2,start )=sci1 703 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 704 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 705 706 case(4) 707 708 scr1=0.0d0 ; sci1=0.0d0 709 scr2=0.0d0 ; sci2=0.0d0 710 scr3=0.0d0 ; sci3=0.0d0 711 scr4=0.0d0 ; sci4=0.0d0 712 do ipw=1,nincpw 713 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 714 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 715 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 716 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 717 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 718 end do 719 scalars(1,start )=scr1 ; scalars(2,start )=sci1 720 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 721 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 722 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 723 724 case(5) 725 726 scr1=0.0d0 ; sci1=0.0d0 727 scr2=0.0d0 ; sci2=0.0d0 728 scr3=0.0d0 ; sci3=0.0d0 729 scr4=0.0d0 ; sci4=0.0d0 730 scr5=0.0d0 ; sci5=0.0d0 731 do ipw=1,nincpw 732 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 733 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 734 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 735 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 736 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 737 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 738 end do 739 scalars(1,start )=scr1 ; scalars(2,start )=sci1 740 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 741 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 742 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 743 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 744 745 case(6) 746 747 scr1=0.0d0 ; sci1=0.0d0 748 scr2=0.0d0 ; sci2=0.0d0 749 scr3=0.0d0 ; sci3=0.0d0 750 scr4=0.0d0 ; sci4=0.0d0 751 scr5=0.0d0 ; sci5=0.0d0 752 scr6=0.0d0 ; sci6=0.0d0 753 do ipw=1,nincpw 754 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 755 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 756 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 757 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 758 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 759 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 760 scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw) 761 end do 762 scalars(1,start )=scr1 ; scalars(2,start )=sci1 763 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 764 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 765 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 766 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 767 scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6 768 769 case(7) 770 771 scr1=0.0d0 ; sci1=0.0d0 772 scr2=0.0d0 ; sci2=0.0d0 773 scr3=0.0d0 ; sci3=0.0d0 774 scr4=0.0d0 ; sci4=0.0d0 775 scr5=0.0d0 ; sci5=0.0d0 776 scr6=0.0d0 ; sci6=0.0d0 777 scr7=0.0d0 ; sci7=0.0d0 778 do ipw=1,nincpw 779 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 780 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 781 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 782 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 783 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 784 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 785 scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw) 786 scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw) 787 end do 788 scalars(1,start )=scr1 ; scalars(2,start )=sci1 789 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 790 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 791 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 792 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 793 scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6 794 scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7 795 796 case(8) 797 798 scr1=0.0d0 ; sci1=0.0d0 799 scr2=0.0d0 ; sci2=0.0d0 800 scr3=0.0d0 ; sci3=0.0d0 801 scr4=0.0d0 ; sci4=0.0d0 802 scr5=0.0d0 ; sci5=0.0d0 803 scr6=0.0d0 ; sci6=0.0d0 804 scr7=0.0d0 ; sci7=0.0d0 805 scr8=0.0d0 ; sci8=0.0d0 806 do ipw=1,nincpw 807 ar=teffv(1,ipw) ; ai=teffv(2,ipw) 808 scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw) 809 scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw) 810 scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw) 811 scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw) 812 scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw) 813 scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw) 814 scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw) 815 scr8=scr8+ar*ffkg(start+7,ipw) ; sci8=sci8+ai*ffkg(start+7,ipw) 816 end do 817 scalars(1,start )=scr1 ; scalars(2,start )=sci1 818 scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2 819 scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3 820 scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4 821 scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5 822 scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6 823 scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7 824 scalars(1,start+7)=scr8 ; scalars(2,start+7)=sci8 825 826 end select 827 828 ! End loop on start 829 end do 830 831 ! End if statement for small or big nffkg 832 end if 833 834 ! ******* Leaving the critical part ********************************* 835 836 ! DEBUG 837 ! write(std_out,*)' opernl4a, write scalars ' 838 ! do iffkg=1,nffkg 839 ! write(std_out,*)iffkg,scalars(1:2,iffkg) 840 ! end do 841 ! ENDDEBUG 842 843 if(istwf_k>=2)then 844 ! Impose parity of resulting scalar (this operation could be 845 ! replaced by direct saving of CPU time in the preceeding section) 846 do iffkg=1,nffkg 847 scalars(parity(iffkg),iffkg)=0.0d0 848 end do 849 end if 850 851 iffkg=0 ; iffkgs=nffkge+nffkgd ; iffkgk=nffkge*2 852 iffkgs2=nffkge+nffkgs 853 do ilang=1,nlang 854 nproj=jproj(ilang) 855 if(nproj>0)then 856 ! ilang2 is the number of independent tensor components 857 ! for symmetric tensor of rank ilang-1 858 ilang2=(ilang*(ilang+1))/2 859 860 ! Loop over projectors 861 do iproj=1,nproj 862 ! Multiply by the k+G factors (tensors of various rank) 863 do ii=1,ilang2 864 ! Get the starting address for the relevant tensor 865 jj=ii+((ilang-1)*ilang*(ilang+1))/6 866 iffkg=iffkg+1 867 ! !$OMP CRITICAL (OPERNL4a_1) 868 gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg) 869 gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg) 870 ! !$OMP END CRITICAL (OPERNL4a_1) 871 ! Now, compute gradients, if needed. 872 if ((choice==2.or.choice==23) .and. ndgxdt==3) then 873 do mu=1,3 874 jffkg=nffkge+(iffkg-1)*3+mu 875 ! Pay attention to the use of reals and imaginary parts here ... 876 ! !$OMP CRITICAL (OPERNL4a_2) 877 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 878 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 879 ! !$OMP END CRITICAL (OPERNL4a_2) 880 end do 881 end if 882 if (choice==2 .and. ndgxdt==1) then 883 jffkg=nffkge+iffkg 884 ! Pay attention to the use of reals and imaginary parts here ... 885 ! !$OMP CRITICAL (OPERNL4a_3) 886 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg) 887 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg) 888 ! !$OMP END CRITICAL (OPERNL4a_3) 889 end if 890 if (choice==4) then 891 do mu=1,3 892 jffkg=nffkge+(iffkg-1)*9+mu 893 ! Pay attention to the use of reals and imaginary parts here ... 894 ! !$OMP CRITICAL (OPERNL4a_4) 895 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 896 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 897 ! !$OMP END CRITICAL (OPERNL4a_4) 898 end do 899 do mu=4,9 900 jffkg=nffkge+(iffkg-1)*9+mu 901 ! Pay attention to the use of reals and imaginary parts here ... 902 ! Also, note the multiplication by (2 pi)**2 903 ! !$OMP CRITICAL (OPERNL4a_5) 904 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg) 905 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg) 906 ! !$OMP END CRITICAL (OPERNL4a_5) 907 end do 908 end if 909 ! End loop on ii=1,ilang2 910 end do 911 912 if ((choice==3.or.choice==23) .or. choice==6) then 913 ! Compute additional tensors related to strain gradients 914 ! ilang3 is number of unique tensor components of rank ilang+1 915 ilang3=((ilang+2)*(ilang+3))/2 916 jjs=((ilang+1)*(ilang+2)*(ilang+3))/6 917 ! Compute strain gradient tensor components 918 do ii=1,ilang3 919 ! Note that iffkgs is also used by ddk and 2nd derivative parts 920 iffkgs=iffkgs+1 921 jj=ii+jjs 922 ! !$OMP CRITICAL (OPERNL4a_6) 923 dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs) 924 dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs) 925 ! !$OMP END CRITICAL (OPERNL4a_6) 926 end do 927 end if 928 929 if (choice==6) then 930 ! Compute additional tensors related to strain 2nd derivatives 931 ! and internal strain derivatives 932 ! ilang6 is number of unique tensor components of rank ilang+3 933 ilang6=((ilang+4)*(ilang+5))/2 934 jjs=((ilang+3)*(ilang+4)*(ilang+5))/6 935 ! Compute strain gradient tensor components 936 do ii=1,ilang6 937 ! Note that iffkgs is also used by ddk part 938 iffkgs2=iffkgs2+1 939 jj=ii+jjs 940 ! !$OMP CRITICAL (OPERNL4a_7) 941 d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2) 942 d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2) 943 ! !$OMP END CRITICAL (OPERNL4a_7) 944 end do 945 946 ! ilang4 is number of unique tensor components of rank ilang 947 ilang4=((ilang+1)*(ilang+2))/2 948 jjs=((ilang)*(ilang+1)*(ilang+2))/6 949 ! Compute internal strain gradient tensor components 950 do ii=1,ilang4 951 iffkgs2=iffkgs2+1 952 jj=ii+jjs 953 ! !$OMP CRITICAL (OPERNL4a_8) 954 ! Pay attention to the use of reals and imaginary parts here ... 955 dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2) 956 dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2) 957 ! !$OMP END CRITICAL (OPERNL4a_8) 958 end do 959 960 ! ilang5 is number of unique tensor components of rank ilang+2 961 ilang5=((ilang+3)*(ilang+4))/2 962 jjs=((ilang+2)*(ilang+3)*(ilang+4))/6 963 ! Compute internal strain gradient tensor components 964 do ii=1,ilang5 965 iffkgs2=iffkgs2+1 966 jj=ii+jjs 967 ! !$OMP CRITICAL (OPERNL4a_9) 968 ! Pay attention to the use of reals and imaginary parts here ... 969 d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2) 970 d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2) 971 ! !$OMP END CRITICAL (OPERNL4a_9) 972 end do 973 end if ! choice==6 974 975 if (choice==5) then 976 ! Compute additional tensors related to ddk with ffnl(:,2,...) 977 ilangx=(ilang*(ilang+1))/2 978 jjs=((ilang-1)*ilang*(ilang+1))/6 979 do ii=1,ilangx 980 ! Note that iffkgs is also used by strain part 981 iffkgs=iffkgs+1 982 jj=ii+jjs 983 ! !$OMP CRITICAL (OPERNL4a_10) 984 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs) 985 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs) 986 ! !$OMP END CRITICAL (OPERNL4a_10) 987 end do 988 ! Compute additional tensors related to ddk with ffnl(:,1,...) 989 if(ilang>=2)then 990 ilangx=((ilang-1)*ilang)/2 991 jjs=((ilang-2)*(ilang-1)*ilang)/6 992 do ii=1,ilangx 993 iffkgk=iffkgk+1 994 jj=ii+jjs 995 ! !$OMP CRITICAL (OPERNL4a_11) 996 dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk) 997 dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk) 998 ! !$OMP END CRITICAL (OPERNL4a_11) 999 end do 1000 end if 1001 end if 1002 1003 ! End projector loop 1004 end do 1005 1006 ! End condition of non-zero projectors 1007 end if 1008 1009 ! End angular momentum loop 1010 end do 1011 1012 ! End loop on atoms 1013 end do 1014 1015 ! End loop on blocks of planewaves 1016 end do 1017 !!$OMP END DO 1018 1019 ABI_DEALLOCATE(ffkg) 1020 ABI_DEALLOCATE(kpgx) 1021 ABI_DEALLOCATE(parity) 1022 ABI_DEALLOCATE(scalars) 1023 ABI_DEALLOCATE(teffv) 1024 !!$OMP END PARALLEL 1025 1026 !DEBUG 1027 !write(std_out,*)' opernl4a : exit' 1028 !ENDDEBUG 1029 1030 end subroutine opernl4a