TABLE OF CONTENTS
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)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
if(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.
PARENTS
nonlop_pl
CHILDREN
dfpt_mkffkg
SOURCE
99 #if defined HAVE_CONFIG_H 100 #include "config.h" 101 #endif 102 103 #include "abi_common.h" 104 105 subroutine opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,& 106 & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,& 107 & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,& 108 & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,& 109 & ntypat,ph3d,sign,vect) 110 111 use defs_basis 112 use m_profiling_abi 113 use m_errors 114 115 !This section has been created automatically by the script Abilint (TD). 116 !Do not modify the following lines by hand. 117 #undef ABI_FUNC 118 #define ABI_FUNC 'opernl3' 119 use interfaces_66_nonlocal, except_this_one => opernl3 120 !End of the abilint section 121 122 implicit none 123 124 !Arguments ------------------------------------ 125 !scalars 126 integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk 127 integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt 128 integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign 129 !arrays 130 integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw) 131 integer,intent(in) :: nloalg(3) 132 real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3) 133 real(dp),intent(in) :: kpg_k(npw,nkpg),kpt(3),ph3d(2,npw,matblk) 134 real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) 135 real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) 136 real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:) 137 real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj) 138 real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj) 139 real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj) 140 141 !Local variables------------------------------- 142 !ntens sets the maximum number of independent tensor components 143 !over all allowed angular momenta; need 20 for spdf for tensors 144 !up to rank 3; to handle stress tensor, need up to rank 5 145 !to handle strain 2DER, need up to rank 7 146 !scalars 147 integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3 148 integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw 149 integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj 150 integer :: ntens 151 real(dp) :: ai,ar 152 real(dp),parameter :: two_pi2=two_pi**2 153 character(len=500) :: message 154 !arrays 155 integer,allocatable :: parity(:) 156 ! real(dp) :: tsec(2) 157 real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:) 158 159 ! ************************************************************************* 160 161 !call wrtout(std_out,"in opernl3","COLL") 162 163 !mblkpw sets the size of blocks of planewaves 164 mblkpw=NLO_MBLKPW 165 166 !Get the actual maximum number of projectors 167 mmproj=maxval(indlmn(3,:,itypat)) 168 169 !Initialisation before blocking on the plane waves 170 171 if (sign==1) then 172 ! Put projected scalars to zero 173 gxa(:,:,:,1:mmproj)=0.0d0 174 if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0 175 if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0 176 if (choice==6) then 177 dgxdis(:,:,:,1:mmproj)=0.0d0 178 d2gxdis(:,:,:,1:mmproj)=0.0d0 179 d2gxds2(:,:,:,1:mmproj)=0.0d0 180 end if 181 end if 182 183 !Set up dimension of kpgx and allocate 184 !ntens sets the maximum number of independent tensor components 185 !over all allowed angular momenta; need 20 for spdf for tensors 186 !up to rank 3; to handle stress tensor, need up to rank 5 187 ntens=1 188 if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4 189 if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10 190 if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20 191 if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35 192 if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56 193 if(choice==6 .and. nlang>=3)ntens=84 194 if(choice==6 .and. nlang==4)ntens=120 195 196 !Set up second dimension of ffkg array, and allocate 197 nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0 198 do ilang=1,nlang 199 ! Get the number of projectors for that angular momentum 200 nproj=jproj(ilang) 201 ! If there is a non-local part, accumulate the number of vectors needed 202 ! The variables ilang below are the number of independent tensors of 203 ! various ranks, the variable names being more historical than logical. 204 ! ilang2=number of rank ilang-1 205 ! ilang3=number of rank ilang+1 206 ! ilang4=number of rank ilang 207 ! ilang5=number of rank ilang+2 208 ! ilang6=number of rank ilang+3 209 if(nproj>0)then 210 ilang2=(ilang*(ilang+1))/2 211 nffkge=nffkge+nproj*ilang2 212 if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang) 213 if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2 214 if(choice==3 .or. choice==6 .or. choice==23)then 215 ilang3=((ilang+2)*(ilang+3))/2 216 nffkgs=nffkgs+nproj*ilang3 217 end if 218 if(choice==6)then 219 ilang4=((ilang+1)*(ilang+2))/2 220 ilang5=((ilang+3)*(ilang+4))/2 221 ilang6=((ilang+4)*(ilang+5))/2 222 nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6) 223 end if 224 end if 225 end do 226 nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk 227 228 !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a 229 !sigfault with nthreads==1 230 ! 231 !Loop on subsets of plane waves (blocking) 232 !!$OMP PARALLEL DEFAULT(PRIVATE) & 233 !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) & 234 !!$OMP SHARED(istwf_k,itypat,jproj,kg_k,kpg_k,kpt,lmnmax,mblkpw,mproj) & 235 !!$OMP SHARED(ndgxdt,nffkg,nffkgd,nffkge,nffkgs,nincat,nkpg,nlang) & 236 !!$OMP SHARED(nloalg,ph3d,npw,ntens,ntypat,sign,vect) 237 238 ABI_ALLOCATE(ffkg,(nffkg,mblkpw)) 239 ABI_ALLOCATE(parity,(nffkg)) 240 ABI_ALLOCATE(kpgx,(mblkpw,ntens)) 241 ABI_ALLOCATE(scalars,(2,nffkg)) 242 ABI_ALLOCATE(teffv,(2,mblkpw)) 243 !!$OMP DO 244 do ipw1=1,npw,mblkpw 245 246 ipw2=min(npw,ipw1+mblkpw-1) 247 nincpw=ipw2-ipw1+1 248 249 ! Initialize kpgx array related to tensors defined below 250 call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,& 251 & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,& 252 & npw,ntens,ntypat,parity) 253 254 ! call timab(74+choice,2,tsec) 255 256 ! Now treat the different signs 257 if (sign==1) then 258 259 do ia=1,nincat 260 261 ! Compute the shift eventually needed to get the phases in ph3d 262 iaph3d=ia 263 if(nloalg(2)>0)iaph3d=ia+ia3-1 264 265 do iffkg=1,nffkg 266 scalars(1,iffkg)=0.0d0 267 scalars(2,iffkg)=0.0d0 268 end do 269 270 ! ******* Entering the first time-consuming part of the routine ******* 271 272 ! Note : first multiply by the phase factor 273 ! This allows to be left with only real operations afterwards. 274 ig=ipw1 275 276 do ipw=1,nincpw 277 ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d) 278 ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d) 279 do iffkg=1,nffkg 280 scalars(1,iffkg)=scalars(1,iffkg)+ar*ffkg(iffkg,ipw) 281 scalars(2,iffkg)=scalars(2,iffkg)+ai*ffkg(iffkg,ipw) 282 end do 283 ig=ig+1 284 end do 285 286 ! ******* Leaving the critical part ********************************* 287 288 ! DEBUG 289 ! write(std_out,*)' opernl2 : scalars' 290 ! do iffkg=1,10 291 ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg) 292 ! end do 293 ! stop 294 ! ENDDEBUG 295 296 if(istwf_k>=2)then 297 ! Impose parity of resulting scalar (this operation could be 298 ! replaced by direct saving of CPU time in the preceeding section) 299 do iffkg=1,nffkg 300 scalars(parity(iffkg),iffkg)=0.0d0 301 end do 302 end if 303 304 iffkg=0 305 iffkgs=nffkge+nffkgd 306 iffkgs2=nffkge+nffkgs 307 iffkgk=nffkge*2 308 do ilang=1,nlang 309 nproj=jproj(ilang) 310 if(nproj>0)then 311 ! ilang2 is the number of independent tensor components 312 ! for symmetric tensor of rank ilang-1 313 ilang2=(ilang*(ilang+1))/2 314 315 ! Loop over projectors 316 do iproj=1,nproj 317 ! Multiply by the k+G factors (tensors of various rank) 318 do ii=1,ilang2 319 ! Get the starting address for the relevant tensor 320 jj=ii+((ilang-1)*ilang*(ilang+1))/6 321 iffkg=iffkg+1 322 ! !$OMP CRITICAL (OPERNL3_1) 323 gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg) 324 gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg) 325 ! !$OMP END CRITICAL (OPERNL3_1) 326 ! Now, compute gradients, if needed. 327 if ((choice==2.or.choice==23) .and. ndgxdt==3) then 328 do mu=1,3 329 jffkg=nffkge+(iffkg-1)*3+mu 330 ! Pay attention to the use of reals and imaginary parts here ... 331 ! !$OMP CRITICAL (OPERNL3_2) 332 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 333 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 334 ! !$OMP END CRITICAL (OPERNL3_2) 335 end do 336 end if 337 if (choice==2 .and. ndgxdt==1) then 338 jffkg=nffkge+iffkg 339 ! Pay attention to the use of reals and imaginary parts here ... 340 ! !$OMP CRITICAL (OPERNL3_3) 341 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg) 342 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg) 343 ! !$OMP END CRITICAL (OPERNL3_3) 344 end if 345 if (choice==4) then 346 do mu=1,3 347 jffkg=nffkge+(iffkg-1)*9+mu 348 ! Pay attention to the use of reals and imaginary parts here ... 349 ! !$OMP CRITICAL (OPERNL3_4) 350 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg) 351 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg) 352 ! !$OMP END CRITICAL (OPERNL3_4) 353 end do 354 do mu=4,9 355 jffkg=nffkge+(iffkg-1)*9+mu 356 ! Pay attention to the use of reals and imaginary parts here ... 357 ! Also, note the multiplication by (2 pi)**2 358 ! !$OMP CRITICAL (OPERNL3_5) 359 dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg) 360 dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg) 361 ! !$OMP END CRITICAL (OPERNL3_5) 362 end do 363 end if 364 ! End loop on ii=1,ilang2 365 end do 366 367 if (choice==3 .or. choice==6 .or. choice==23) then 368 ! Compute additional tensors related to strain gradients 369 ! ilang3 is number of unique tensor components of rank ilang+1 370 ilang3=((ilang+2)*(ilang+3))/2 371 jjs=((ilang+1)*(ilang+2)*(ilang+3))/6 372 ! Compute strain gradient tensor components 373 do ii=1,ilang3 374 ! Note that iffkgs is also used by ddk and 2nd derivative parts 375 iffkgs=iffkgs+1 376 jj=ii+jjs 377 ! !$OMP CRITICAL (OPERNL3_6) 378 dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs) 379 dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs) 380 ! !$OMP END CRITICAL (OPERNL3_6) 381 end do 382 end if 383 384 if (choice==6) then 385 ! Compute additional tensors related to strain 2nd derivatives 386 ! and internal strain derivatives 387 ! ilang6 is number of unique tensor components of rank ilang+3 388 ilang6=((ilang+4)*(ilang+5))/2 389 jjs=((ilang+3)*(ilang+4)*(ilang+5))/6 390 ! Compute strain gradient tensor components 391 do ii=1,ilang6 392 ! Note that iffkgs is also used by ddk part 393 iffkgs2=iffkgs2+1 394 jj=ii+jjs 395 ! !$OMP CRITICAL (OPERNL3_6) 396 d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2) 397 d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2) 398 ! !$OMP END CRITICAL (OPERNL3_6) 399 end do 400 401 ! ilang4 is number of unique tensor components of rank ilang 402 ilang4=((ilang+1)*(ilang+2))/2 403 jjs=((ilang)*(ilang+1)*(ilang+2))/6 404 ! Compute internal strain gradient tensor components 405 do ii=1,ilang4 406 iffkgs2=iffkgs2+1 407 jj=ii+jjs 408 ! !$OMP CRITICAL (OPERNL3_6) 409 ! Pay attention to the use of reals and imaginary parts here ... 410 dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2) 411 dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2) 412 ! !$OMP END CRITICAL (OPERNL3_6) 413 end do 414 415 ! ilang5 is number of unique tensor components of rank ilang+2 416 ilang5=((ilang+3)*(ilang+4))/2 417 jjs=((ilang+2)*(ilang+3)*(ilang+4))/6 418 ! Compute internal strain gradient tensor components 419 do ii=1,ilang5 420 iffkgs2=iffkgs2+1 421 jj=ii+jjs 422 ! !$OMP CRITICAL (OPERNL3_6) 423 ! Pay attention to the use of reals and imaginary parts here ... 424 d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2) 425 d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2) 426 ! !$OMP END CRITICAL (OPERNL3_6) 427 end do 428 end if ! choice==6 429 430 if (choice==5) then 431 ! Compute additional tensors related to ddk with ffnl(:,2,...) 432 ilangx=(ilang*(ilang+1))/2 433 jjs=((ilang-1)*ilang*(ilang+1))/6 434 do ii=1,ilangx 435 ! Note that iffkgs is also used by strain part 436 iffkgs=iffkgs+1 437 jj=ii+jjs 438 ! !$OMP CRITICAL (OPERNL3_7) 439 dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs) 440 dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs) 441 ! !$OMP END CRITICAL (OPERNL3_7) 442 end do 443 ! Compute additional tensors related to ddk with ffnl(:,1,...) 444 if(ilang>=2)then 445 ilangx=((ilang-1)*ilang)/2 446 jjs=((ilang-2)*(ilang-1)*ilang)/6 447 do ii=1,ilangx 448 iffkgk=iffkgk+1 449 jj=ii+jjs 450 ! !$OMP CRITICAL (OPERNL3_8) 451 dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk) 452 dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk) 453 ! !$OMP END CRITICAL (OPERNL3_8) 454 end do 455 end if 456 end if 457 458 ! End projector loop 459 end do 460 461 ! End condition of non-zero projectors 462 end if 463 464 ! End angular momentum loop 465 end do 466 467 ! End loop on atoms 468 end do 469 470 else if (sign==-1) then 471 ! Application of non-local part from projected scalars 472 ! back to reciprocal space ... 473 ! [this section merely computes terms which add to <G|Vnl|C>; 474 ! nothing here is needed when various gradients are being computed] 475 476 ! Loop on atoms 477 do ia=1,nincat 478 479 ! Compute the shift eventually needed to get the phases in ph3d 480 iaph3d=ia 481 if(nloalg(2)>0)iaph3d=ia+ia3-1 482 483 ! Transfer gxa (and eventually dgxdt) in scalars with different indexing 484 iffkg=0 485 iffkgk=nffkge*2 486 iffkgs=nffkge 487 do ilang=1,nlang 488 nproj=jproj(ilang) 489 if (nproj>0) then 490 ilang2=(ilang*(ilang+1))/2 491 ilang3=((ilang+2)*(ilang+3))/2 492 do iproj=1,nproj 493 do ii=1,ilang2 494 jj=ii+((ilang-1)*ilang*(ilang+1))/6 495 iffkg=iffkg+1 496 if(choice==1 .or. choice==3)then 497 scalars(1,iffkg)=gxa(1,jj,ia,iproj) 498 scalars(2,iffkg)=gxa(2,jj,ia,iproj) 499 else if (choice==2 .and. ndgxdt==1) then 500 jffkg=nffkge+iffkg 501 ! Pay attention to the use of reals and imaginary parts here ... 502 ! Also, the gxa and dgxdt arrays are switched, in order 503 ! to give the correct combination when multiplying ffkg, 504 ! see Eq.(53) of PRB55,10337(1997) 505 scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj) 506 scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj) 507 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 508 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 509 else if (choice==5) then 510 jffkg=nffkge+iffkg 511 ! The gxa and dgxdt arrays are switched, in order 512 ! to give the correct combination when multiplying ffkg, 513 scalars(1,jffkg)= gxa(1,jj,ia,iproj) 514 scalars(2,jffkg)= gxa(2,jj,ia,iproj) 515 scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj) 516 scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj) 517 end if 518 end do 519 if(choice==3) then 520 do ii=1,ilang3 521 iffkgs=iffkgs+1 522 jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6 523 scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj) 524 scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj) 525 end do 526 end if 527 if(ilang>=2 .and. choice==5)then 528 do ii=1,((ilang-1)*ilang)/2 529 jj=ii+((ilang-2)*(ilang-1)*ilang)/6 530 iffkgk=iffkgk+1 531 scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj) 532 scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj) 533 end do 534 end if 535 end do 536 end if 537 end do 538 539 ! DEBUG 540 ! if(choice==5)then 541 ! write(std_out,*)' opernl2 : write gxa(:,...) array ' 542 ! do ii=1,10 543 ! write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1) 544 ! end do 545 ! write(std_out,*)' opernl2 : write dgxdt(:,1,...) array ' 546 ! do ii=1,10 547 ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1) 548 ! end do 549 ! write(std_out,*)' opernl2 : write dgxdt(:,2,...) array ' 550 ! do ii=1,4 551 ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1) 552 ! end do 553 ! end if 554 555 ! do iffkg=1,nffkg 556 ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg) 557 ! end do 558 ! stop 559 ! ENDDEBUG 560 561 ig=ipw1 562 563 ! ******* Entering the critical part ********************************* 564 565 do ipw=1,nincpw 566 ar=0.0d0 567 ai=0.0d0 568 do iffkg=1,nffkg 569 ar=ar+ffkg(iffkg,ipw)*scalars(1,iffkg) 570 ai=ai+ffkg(iffkg,ipw)*scalars(2,iffkg) 571 end do 572 vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d) 573 vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d) 574 ig=ig+1 575 end do 576 577 ! ******* Leaving the critical part ********************************* 578 579 ! End loop on atoms 580 end do 581 582 ! End sign==-1 583 else 584 585 ! Problem: sign and/or choice do not make sense 586 write(message, '(a,2i10,a,a)' )& 587 & ' Input sign, choice=',sign,choice,ch10,& 588 & ' Are not compatible or allowed. ' 589 MSG_BUG(message) 590 end if 591 592 ! End loop on blocks of planewaves 593 end do 594 !!$OMP END DO 595 ABI_DEALLOCATE(ffkg) 596 ABI_DEALLOCATE(kpgx) 597 ABI_DEALLOCATE(parity) 598 ABI_DEALLOCATE(scalars) 599 ABI_DEALLOCATE(teffv) 600 !!$OMP END PARALLEL 601 602 603 !DEBUG 604 !if(choice==5)then 605 !write(std_out,*)'opernl2 : write vect(2*npw)' 606 !do ipw=1,2 607 !write(std_out,*)ipw,vect(1:2,ipw) 608 !end do 609 !write(std_out,*)'opernl2 : write ph3d' 610 !do ipw=1,npw 611 !write(std_out,*)ipw,ph3d(1:2,ipw,1) 612 !end do 613 !write(std_out,*)' opernl2 : write gxa array ' 614 !write(std_out,*)' ang mom ,ia ' 615 !do iproj=1,mproj 616 !do ia=1,1 617 !do ii=1,3 618 !do ii=1,10 619 !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1) 620 !end do 621 !end do 622 !end do 623 !end if 624 !if(choice==5)then 625 !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array ' 626 !do ii=1,10 627 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1) 628 !end do 629 !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array ' 630 !do ii=1,4 631 !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1) 632 !end do 633 !stop 634 !end if 635 !ENDDEBUG 636 637 end subroutine opernl3