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