TABLE OF CONTENTS
ABINIT/m_opernla_ylm [ Modules ]
NAME
m_opernla_ylm
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MT) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_opernla_ylm 22 23 use defs_basis 24 use m_abicore 25 use m_errors 26 use m_xmpi 27 #if defined HAVE_OPENMP 28 use OMP_LIB 29 #endif 30 31 use defs_abitypes, only : MPI_type 32 use m_time, only : timab 33 34 implicit none 35 36 private
ABINIT/opernla_ylm [ Functions ]
NAME
opernla_ylm
FUNCTION
For a given wave-function |c>, get all projected scalars <p_lmn|c> where |p_lmn> are non-local projectors With: <p_lmn|c>=4pi/sqrt(vol) (i)^l Sum_g[c(g).f_nl(g).Y_lm(g).exp(2pi.i.g.R)]
INPUTS
choice=chooses possible output:
if choice>=0: compute projected scalars
if choice<0: same as choice>0 but use already computed projected scalars
if ABS(choice)>1, then compute additional quantities:
2: compute projected scalars and derivatives wrt atm pos.
3: compute projected scalars and derivatives wrt strains
22: compute projected scalars and 2nd derivatives wrt atm pos. and q-vector.
23: compute projected scalars, derivatives wrt atm pos. and derivatives wrt strains
25: compute projected scalars and 3rd derivatives wrt atm pos. and two q-vectors.
4, 24: compute projected scalars, derivatives wrt atm pos.
and 2nd derivatives wrt atm pos.
33: compute projected scalars and 2nd derivatives wrt strain and q-vector.
5,51,52: compute projected scalars and derivatives wrt wave vector k
53: compute projected scalars and derivatives wrt wave vector k in direction idir+1 and idir+2 mod 3
54: compute projected scalars, deriv. wrt atm pos., deriv. wrt wave vector k
and 2nd derivatives wrt right wave vector k and atm pos.
55: compute projected scalars, deriv. strains, deriv. wrt wave vector k
and 2nd derivatives wrt right wave vector k and strain
6: compute projected scalars, derivatives wrt atm pos., derivatives wrt strains,
2nd derivatives wrt 2 strains and derivatives wrt strain and atm pos.
7: not available
8: compute projected scalars, derivatives wrt wave vector k
and 2nd derivatives wrt 2 wave vectors k
cplex=1 if <p_lmn|c> scalars are real or pure imaginary (equivalent to istwfk>1)
2 if <p_lmn|c> scalars are complex
dimffnl=second dimension of ffnl
ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
ia3=gives the number of the first atom in the subset presently treated
idir=direction of the - atom to be moved in the case (choice=2,signs=2) or (choice=22,signs=2)
- k point direction in the case (choice=5,signs=2)
- strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
- strain component (1:9) in the case (choice=33,signs=2)
- (1:9) components to specify the atom to be moved and the second q-gradient
direction in the case (choice=25,signs=2)
indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
istwf_k=option parameter that describes the storage of wfs
kpg(npw,nkpg)=(k+G) components for ikpg=1...3 (if nkpg=3 or 9)
[(k+G)_a].[(k+G)_b] quantities for ikpg=4...9 (if nkpg=9)
(k+G) Cartesian components for choice==33
matblk=dimension of the array ph3d
mpi_enreg=information about MPI parallelization
ndgxdt=second dimension of dgxdt
nd2gxdt=second dimension of d2gxdt
nincat=number of atoms in the subset here treated
nkpg=second dimension of array kpg (0, 3 or 9)
nlmn=number of (l,m,n) numbers for current type of atom
nloalg(3)=governs the choice of the algorithm for non-local operator.
npw=number of plane waves in reciprocal space
nspinor=number of spinorial components of the wavefunctions (on current proc)
ph3d(2,npw,matblk)=three-dimensional phase factors
[qdir]= optional, direction of the q-gradient (only for choice=22 choice=25 and choice=33)
signs=chooses possible output:
signs=1: compute derivatives in all directions
signs=2: compute derivative in direction IDIR only
compatible only with 1st-order derivatives and "single" derivatives
ucvol=unit cell volume (bohr^3)
vect(2,npw*my_nspinor)=starting vector in reciprocal space
OUTPUT
if (choice>1) dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)=
gradients of projected scalars wrt coords (choice=2, 23, 4, 54, 6)
wrt strains (choice=3, 23, 55)
wrt k wave vect. (choice=5, 51, 52, 53, 54, 55, 8)
wrt coords and q vect (choice=22)
wrt coords and two q vects (choice=25)
wrt strains and q vect (choice=33)
if (choice=4, 24, 33, 54, 55, 6, 8) d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)=
2nd grads of projected scalars wrt 2 coords (choice=4 or 24)
wrt coords & k wave vect. (choice=54)
wrt strains & k wave vect. (choice=55)
wrt coords & strains (choice=6)
wrt 2 strains (choice=6)
wrt 2 k wave vect. (choice=8)
wrt strains and q vect (choice=33)
only compatible with signs=1
cplex_dgxdt(ndgxdt) = used only when cplex = 1
cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
cplex_d2gxdt(nd2gxdt) = used only when cplex = 1
cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary
SIDE EFFECTS
gx(cplex,nlmn,nincat,nspinor)= projected scalars - input if choice<0, output if choice>=0
NOTES
1-The openMP version is different from the standard version: the standard version is more effifient on one CPU core. 2-Operate for one type of atom, and within this given type of atom, for a subset of at most nincat atoms.
SOURCE
147 subroutine opernla_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,d2gxdt,dgxdt,ffnl,gx,& 148 & ia3,idir,indlmn,istwf_k,kpg,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpg,nlmn,& 149 & nloalg,npw,nspinor,ph3d,signs,ucvol,vect,qdir) 150 151 !Arguments ------------------------------------ 152 !scalars 153 integer,intent(in) :: choice,cplex,dimffnl,ia3,idir,istwf_k,matblk,nd2gxdt 154 integer,intent(in) :: ndgxdt,nincat,nkpg,nlmn,npw,nspinor,signs 155 integer,intent(in),optional :: qdir 156 real(dp),intent(in) :: ucvol 157 type(MPI_type),intent(in) :: mpi_enreg 158 !arrays 159 integer,intent(in) :: indlmn(6,nlmn),nloalg(3) 160 integer,intent(out) :: cplex_dgxdt(ndgxdt),cplex_d2gxdt(nd2gxdt) 161 real(dp),intent(in) :: ffnl(npw,dimffnl,nlmn),kpg(npw,nkpg),ph3d(2,npw,matblk) 162 real(dp),intent(in) :: vect(:,:) 163 real(dp),intent(out) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor) 164 real(dp),intent(out) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor),gx(cplex,nlmn,nincat,nspinor) 165 166 !Local variables------------------------------- 167 !scalars 168 integer :: choice_,gama,gamb,gamc,gamd,ia,iaph3d,ialpha,ibeta,idelta,idelgam 169 integer :: ierr,igamma,il,ilmn,ipw,ipw0,ipwshft,ispinor,jpw,mu,mu0 170 integer :: mua,mub,ndir,nthreads,nua1,nua2,nub1,nub2 171 real(dp), parameter :: two_pi2=two_pi*two_pi 172 real(dp) :: aux_i,aux_i2,aux_i3,aux_i4 173 real(dp) :: aux_r,aux_r2,aux_r3,aux_r4 174 real(dp) :: buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5,buffer_i6 175 real(dp) :: buffer_ia,buffer_ib,buffer_ic,buffer_id,buffer_ie,buffer_if 176 real(dp) :: buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5,buffer_r6 177 real(dp) :: buffer_ra,buffer_rb,buffer_rc,buffer_rd,buffer_re,buffer_rf 178 real(dp) :: kpga,kpgb,kpgc,kpgd,scale,scale2,wt 179 logical :: check,parity 180 !arrays 181 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 182 integer,parameter :: nalpha(9)=(/1,2,3,3,3,2,2,1,1/),nbeta(9)=(/1,2,3,2,1,1,3,3,2/) 183 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 184 integer,parameter :: ffnl_dir_dat(6)=(/2,3,3,1,1,2/) 185 integer,parameter :: idir1(9)=(/1,1,1,2,2,2,3,3,3/),idir2(9)=(/1,2,3,1,2,3,1,2,3/) 186 integer :: ffnl_dir(2) 187 real(dp) :: tsec(2) 188 ! real(dp),allocatable :: kpg(:,:) 189 real(dp),allocatable :: scali(:),scalr(:),scalari(:,:),scalarr(:,:) 190 ! ************************************************************************* 191 192 if (choice==-1) return 193 194 !Useful variables 195 choice_=abs(choice) 196 wt=four_pi/sqrt(ucvol);if (cplex==1) wt=2.d0*wt 197 ipw0=1;if (istwf_k==2.and.mpi_enreg%me_g0_fft==1) ipw0=2 198 cplex_dgxdt(:) = 0 ; if (cplex == 1) cplex_dgxdt(:) = 1 199 cplex_d2gxdt(:) = 0 ; if (cplex == 1) cplex_d2gxdt(:) = 1 200 nthreads=1 201 202 !Check compatibility of options 203 if (signs==1) then 204 ! signs=1, almost all choices 205 check=(choice_==0 .or.choice_==1 .or.choice_==2 .or.choice_==3 .or.choice_==4 .or.& 206 & choice_==23.or.choice_==24.or.choice_==5 .or.choice_==51.or.choice_==52.or.& 207 & choice_==53.or.choice_==54.or.choice_==55.or.& 208 & choice_==6 .or.choice_==8 .or.choice_==81) 209 ABI_CHECK(check,'BUG: choice not compatible (for signs=1)') 210 end if 211 if (signs==2) then 212 ! signs=2,less choices 213 check=(choice_== 0.or.choice_== 1.or.choice_== 2.or.choice_==3.or.choice_==5.or.& 214 & choice_==22.or.choice_==23.or.choice==25.or.choice_==33.or.& 215 & choice_==51.or.choice_==52.or.choice_==53.or.choice_==54.or.choice_==55.or.& 216 & choice_== 8.or.choice_==81) 217 ABI_CHECK(check,'BUG: signs=2 not compatible with this choice') 218 end if 219 if (choice_==3.and.signs==2) then 220 ! 1<=idir<=6 is required when choice=3 and signs=2 221 check=(idir>=1.and.idir<=6) 222 ABI_CHECK(check,'BUG: choice=3 and signs=2 requires 1<=idir<=6') 223 else if (choice_==25.and.signs==2) then 224 ! 1<=idir<=9 is required when choice=25 and signs=2 225 check=(idir>=1.and.idir<=9) 226 ABI_CHECK(check,'BUG: choice=25 and signs=2 requires 1<=idir<=9') 227 else if (choice_==33.and.signs==2) then 228 ! 1<=idir<=9 is required when choice=33 and signs=2 229 check=(idir>=1.and.idir<=9) 230 ABI_CHECK(check,'BUG: choice=33 and signs=2 requires 1<=idir<=9') 231 else if ((choice_==8.or.choice_==81.or.choice_==54).and.signs==2) then 232 ! 1<=idir<=9 is required when choice=8 and signs=2 233 check=(idir>=1.and.idir<=9) 234 ABI_CHECK(check,'BUG: choice=8/81 and signs=2 requires 1<=idir<=9') 235 else 236 ! signs=2 requires 1<=idir<=3 when choice>1 237 check=(signs/=2.or.choice<=1.or.(idir>=1.and.idir<=3)) 238 ABI_CHECK(check,'BUG: signs=2 requires 1<=idir<=3') 239 end if 240 241 #if defined HAVE_OPENMP 242 nthreads=OMP_GET_NUM_THREADS() 243 #endif 244 245 if (opernla_counter>=0) then 246 if (choice>=0) opernla_counter = opernla_counter + 1 247 if ((signs==1).and.(choice_==2.or.choice_==4.or.choice_==6.or.choice_==23.or.choice_==24.or.choice_==54)) then 248 opernla_counter = opernla_counter + 3 ! 3 directions 249 end if 250 if ((signs==2).and.(choice_==2)) opernla_counter = opernla_counter + 1 ! 1 direction 251 if ((signs==1).and.(choice_==3.or.choice_==6.or.choice_==23.or.choice_==55)) opernla_counter = opernla_counter + 6 ! 6 directions 252 if ((signs==2).and.(choice_==3)) opernla_counter = opernla_counter + 1 ! 1 direction 253 end if 254 255 !========================================================================== 256 !========== STANDARD VERSION ============================================== 257 !========================================================================== 258 !Parallelization OpenMP on nlmn loops 259 if (nthreads==1) then 260 261 ! Loop on spinorial components 262 do ispinor =1,nspinor 263 ipwshft=(ispinor-1)*npw 264 265 ! Allocate work space 266 ABI_MALLOC(scali,(npw)) 267 ABI_MALLOC(scalr,(npw)) 268 269 ! Loop on atoms (blocking) 270 do ia=1,nincat 271 iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1 272 ! Compute c(g).exp(2pi.i.g.R) 273 do ipw=ipw0,npw 274 jpw=ipw+ipwshft 275 scalr(ipw)=(vect(1,jpw)*ph3d(1,ipw,iaph3d)-vect(2,jpw)*ph3d(2,ipw,iaph3d)) 276 scali(ipw)=(vect(2,jpw)*ph3d(1,ipw,iaph3d)+vect(1,jpw)*ph3d(2,ipw,iaph3d)) 277 end do 278 279 if (ipw0==2) then 280 scalr(1)=half*vect(1,1+ipwshft)*ph3d(1,1,iaph3d) 281 scali(1)=half*vect(1,1+ipwshft)*ph3d(2,1,iaph3d) 282 end if 283 284 ! -------------------------------------------------------------------- 285 ! ALL CHOICES: 286 ! Accumulate Gx 287 ! -------------------------------------------------------------------- 288 289 if (choice>=0) then 290 291 do ilmn=1,nlmn 292 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 293 scale=wt;if (il>1) scale=-scale 294 if (cplex==2) then 295 buffer_r1 = zero ; buffer_i1 = zero 296 do ipw=1,npw 297 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn) 298 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn) 299 end do 300 if (parity) then 301 gx(1,ilmn,ia,ispinor) = scale*buffer_r1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_i1 302 else 303 gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_r1 304 end if 305 else 306 if (parity) then 307 buffer_r1 = zero 308 do ipw=1,npw 309 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn) 310 end do 311 gx(1,ilmn,ia,ispinor) = scale*buffer_r1 312 else 313 buffer_i1 = zero 314 do ipw=1,npw 315 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn) 316 end do 317 gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 318 end if 319 end if 320 end do 321 end if 322 323 ! -------------------------------------------------------------------- 324 ! CHOICE= 2, 4, 6, 23, 24, 54 -- SIGNS= 1 325 ! Accumulate dGxdt --- derivative wrt atm pos. --- for all directions 326 ! -------------------------------------------------------------------- 327 if ((signs==1).and. & 328 & (choice_==2.or.choice_==4.or.choice_==6.or.choice_==23.or.choice_==24.or.choice_==54)) then 329 mu0=1;if (choice_==23.or.choice_==6) mu0=7 330 331 do ilmn=1,nlmn 332 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 333 scale=wt;if (il>1) scale=-scale 334 scale2=two_pi*scale 335 if (cplex==2) then 336 buffer_r1 = zero ; buffer_r2 = zero 337 buffer_r3 = zero ; buffer_i1 = zero 338 buffer_i2 = zero ; buffer_i3 = zero 339 do ipw=1,npw 340 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 341 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 342 buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1) 343 buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2) 344 buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3) 345 buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1) 346 buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2) 347 buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3) 348 end do 349 if (parity) then 350 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1 351 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_i1 352 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2 353 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_i2 354 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3 355 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_i3 356 else 357 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_i1 358 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1 359 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2 360 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2 361 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3 362 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3 363 end if 364 else 365 if (parity) then 366 buffer_r1 = zero ; buffer_r2 = zero ; buffer_r3 = zero 367 do ipw=1,npw 368 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 369 buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1) 370 buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2) 371 buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3) 372 end do 373 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1 374 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2 375 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3 376 else 377 buffer_i1 = zero ; buffer_i2 = zero ; buffer_i3 = zero 378 do ipw=1,npw 379 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 380 buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1) 381 buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2) 382 buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3) 383 end do 384 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_i1 385 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2 386 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3 387 end if 388 end if 389 end do 390 end if 391 392 ! -------------------------------------------------------------------- 393 ! CHOICE= 2, -- SIGNS= 2 394 ! Accumulate dGxdt --- derivative wrt atm pos. --- for direction IDIR 395 ! -------------------------------------------------------------------- 396 if ((signs==2).and.(choice_==2)) then 397 398 mu0=1 399 do ilmn=1,nlmn 400 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 401 scale=wt;if (il>1) scale=-scale 402 scale2=two_pi*scale 403 if (cplex==2) then 404 buffer_r1 = zero ; buffer_i1 = zero 405 do ipw=1,npw 406 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 407 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 408 buffer_r1 = buffer_r1 - aux_i*kpg(ipw,idir) 409 buffer_i1 = buffer_i1 + aux_r*kpg(ipw,idir) 410 end do 411 if (parity) then 412 dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 413 dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_i1 414 else 415 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 416 dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 417 end if 418 else 419 if (parity) then 420 buffer_r1 = zero 421 do ipw=1,npw 422 buffer_r1 = buffer_r1 - scali(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir) 423 end do 424 dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 425 else 426 buffer_i1 = zero 427 do ipw=1,npw 428 buffer_i1 = buffer_i1 + scalr(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir) 429 end do 430 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 431 end if 432 end if 433 end do 434 end if 435 436 ! -------------------------------------------------------------------- 437 ! CHOICE= 22 and 25 -- SIGNS= 2 438 ! Accumulate dGxdt --- derivative wrt atm pos. and one or two q-vectors--- 439 ! for directions fixed by IDIR and qdir 440 ! -------------------------------------------------------------------- 441 if (signs==2.and.(choice_==22.or.choice_==25)) then 442 !Use same notation as the notes for clarity 443 if(choice_==22) ialpha=idir 444 if(choice_==25) ialpha=nalpha(idir) 445 446 mu0=1 447 do ilmn=1,nlmn 448 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 449 scale=wt;if (il>1) scale=-scale 450 if (cplex==2) then 451 buffer_r1 = zero ; buffer_i1 = zero 452 do ipw=1,npw 453 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 454 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 455 buffer_r1 = buffer_r1 + aux_r*kpg(ipw,ialpha) 456 buffer_i1 = buffer_i1 + aux_i*kpg(ipw,ialpha) 457 end do 458 if (parity) then 459 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 460 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1 461 else 462 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 463 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1 464 end if 465 else 466 if (parity) then 467 buffer_r1 = zero 468 do ipw=1,npw 469 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,ialpha) 470 end do 471 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 472 else 473 buffer_i1 = zero 474 do ipw=1,npw 475 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,ialpha) 476 end do 477 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 478 end if 479 end if 480 end do 481 end if 482 483 ! -------------------------------------------------------------------- 484 ! CHOICE= 3, 6, 23, 55 -- SIGNS= 1 485 ! Accumulate dGxdt --- derivative wrt strain --- for all directions 486 ! -------------------------------------------------------------------- 487 if ((signs==1).and.(choice_==3.or.choice_==6.or.choice_==23.or.choice_==55)) then 488 489 mu0=1 490 do ilmn=1,nlmn 491 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 492 scale=wt;if (il>1) scale=-scale 493 scale2=half*scale 494 if (cplex==2) then 495 buffer_r1 = zero ; buffer_r2 = zero 496 buffer_r3 = zero ; buffer_r4 = zero 497 buffer_r5 = zero ; buffer_r6 = zero 498 buffer_i1 = zero ; buffer_i2 = zero 499 buffer_i3 = zero ; buffer_i4 = zero 500 buffer_i5 = zero ; buffer_i6 = zero 501 do ipw=1,npw 502 aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn) 503 aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn) 504 aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn) 505 aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn) 506 aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn) 507 aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn) 508 buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1) 509 buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2) 510 buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3) 511 buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3) 512 buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3) 513 buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2) 514 buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1) 515 buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2) 516 buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3) 517 buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3) 518 buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3) 519 buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2) 520 end do 521 if (parity) then 522 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1 523 dgxdt(2,mu0 ,ilmn,ia,ispinor) =-scale*buffer_i1 524 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2 525 dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2 526 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3 527 dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3 528 dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4 529 dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_i4 530 dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5 531 dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_i5 532 dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6 533 dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_i6 534 else 535 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 536 dgxdt(2,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1 537 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 538 dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2 539 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 540 dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3 541 dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4 542 dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4 543 dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5 544 dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5 545 dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6 546 dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6 547 end if 548 else 549 if (parity) then 550 buffer_r1 = zero ; buffer_r2 = zero 551 buffer_r3 = zero ; buffer_r4 = zero 552 buffer_r5 = zero ; buffer_r6 = zero 553 do ipw=1,npw 554 aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn) 555 aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn) 556 aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn) 557 buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1) 558 buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2) 559 buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3) 560 buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3) 561 buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3) 562 buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2) 563 end do 564 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1 565 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2 566 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3 567 dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4 568 dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5 569 dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6 570 else 571 buffer_i1 = zero ; buffer_i2 = zero 572 buffer_i3 = zero ; buffer_i4 = zero 573 buffer_i5 = zero ; buffer_i6 = zero 574 do ipw=1,npw 575 aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn) 576 aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn) 577 aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn) 578 buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1) 579 buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2) 580 buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3) 581 buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3) 582 buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3) 583 buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2) 584 end do 585 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 586 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 587 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 588 dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4 589 dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5 590 dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6 591 end if 592 end if 593 end do 594 end if 595 596 ! -------------------------------------------------------------------- 597 ! CHOICE= 3 -- SIGNS= 2 598 ! Accumulate dGxdt --- derivative wrt strain --- for direction IDIR 599 ! -------------------------------------------------------------------- 600 if ((signs==2).and.(choice_==3)) then 601 602 mu0=1;ffnl_dir(1)=1;if (dimffnl>2) ffnl_dir(1)=idir 603 do ilmn=1,nlmn 604 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 605 scale=wt;if (il>1) scale=-scale 606 if (cplex==2) then 607 buffer_r1 = zero ; buffer_i1 = zero 608 do ipw=1,npw 609 buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 610 buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 611 end do 612 if (parity) then 613 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 614 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1 615 else 616 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 617 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1 618 end if 619 else 620 if (parity) then 621 buffer_r1 = zero 622 do ipw=1,npw 623 buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 624 end do 625 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 626 else 627 buffer_i1 = zero 628 do ipw=1,npw 629 buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 630 end do 631 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 632 end if 633 end if 634 end do 635 end if 636 637 ! -------------------------------------------------------------------- 638 ! CHOICE= 33 -- SIGNS= 2 639 ! Accumulate dGxdt and d2Gxdt --- mixed derivative wrt strain and q-vector --- 640 ! for directions fixed by idir and qdir 641 ! -------------------------------------------------------------------- 642 if ((signs==2).and.(choice_==33)) then 643 644 !Use same notation as the notes for clarity 645 ibeta=nalpha(idir) 646 idelta=nbeta(idir) 647 igamma=qdir 648 idelgam=gamma(idelta,igamma) 649 do ilmn=1,nlmn 650 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 651 scale=wt;if (il>1) scale=-scale 652 if (cplex==2) then 653 buffer_r1 = zero ; buffer_i1 = zero 654 buffer_r2 = zero ; buffer_i2 = zero 655 buffer_r3 = zero ; buffer_i3 = zero 656 buffer_r4 = zero ; buffer_i4 = zero 657 buffer_r5 = zero ; buffer_i5 = zero 658 do ipw=1,npw 659 aux_r = scalr(ipw)*ffnl(ipw,1+igamma,ilmn) 660 aux_i = scali(ipw)*ffnl(ipw,1+igamma,ilmn) 661 buffer_r1 = buffer_r1 + aux_r 662 buffer_i1 = buffer_i1 + aux_i 663 buffer_r2 = buffer_r2 + aux_r*kpg(ipw,ibeta) 664 buffer_i2 = buffer_i2 + aux_i*kpg(ipw,ibeta) 665 aux_r = scalr(ipw)*ffnl(ipw,1+idelta,ilmn) 666 aux_i = scali(ipw)*ffnl(ipw,1+idelta,ilmn) 667 buffer_r3 = buffer_r3 + aux_r 668 buffer_i3 = buffer_i3 + aux_i 669 buffer_r4 = buffer_r4 + aux_r*kpg(ipw,ibeta) 670 buffer_i4 = buffer_i4 + aux_i*kpg(ipw,ibeta) 671 aux_r = scalr(ipw)*ffnl(ipw,4+idelgam,ilmn) 672 aux_i = scali(ipw)*ffnl(ipw,4+idelgam,ilmn) 673 buffer_r5 = buffer_r5 + aux_r*kpg(ipw,ibeta) 674 buffer_i5 = buffer_i5 + aux_i*kpg(ipw,ibeta) 675 end do 676 677 if (parity) then 678 dgxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r1 679 dgxdt(2,1,ilmn,ia,ispinor) = scale*buffer_i1 680 dgxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r3 681 dgxdt(2,2,ilmn,ia,ispinor) = scale*buffer_i3 682 d2gxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r2 683 d2gxdt(2,1,ilmn,ia,ispinor) = scale*buffer_i2 684 d2gxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r4 685 d2gxdt(2,2,ilmn,ia,ispinor) = scale*buffer_i4 686 d2gxdt(1,3,ilmn,ia,ispinor) = scale*buffer_r5 687 d2gxdt(2,3,ilmn,ia,ispinor) = scale*buffer_i5 688 else 689 dgxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i1 690 dgxdt(2,1,ilmn,ia,ispinor) = scale*buffer_r1 691 dgxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i3 692 dgxdt(2,2,ilmn,ia,ispinor) = scale*buffer_r3 693 d2gxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i2 694 d2gxdt(2,1,ilmn,ia,ispinor) = scale*buffer_r2 695 d2gxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i4 696 d2gxdt(2,2,ilmn,ia,ispinor) = scale*buffer_r4 697 d2gxdt(1,3,ilmn,ia,ispinor) =-scale*buffer_i5 698 d2gxdt(2,3,ilmn,ia,ispinor) = scale*buffer_r5 699 end if 700 else 701 if (parity) then 702 buffer_r1 = zero 703 buffer_r2 = zero 704 buffer_r3 = zero 705 buffer_r4 = zero 706 buffer_r5 = zero 707 do ipw=1,npw 708 aux_r = scalr(ipw)*ffnl(ipw,1+igamma,ilmn) 709 buffer_r1 = buffer_r1 + aux_r 710 buffer_r2 = buffer_r2 + aux_r*kpg(ipw,ibeta) 711 aux_r = scalr(ipw)*ffnl(ipw,1+idelta,ilmn) 712 buffer_r3 = buffer_r3 + aux_r 713 buffer_r4 = buffer_r4 + aux_r*kpg(ipw,ibeta) 714 aux_r = scalr(ipw)*ffnl(ipw,4+idelgam,ilmn) 715 buffer_r5 = buffer_r5 + aux_r*kpg(ipw,ibeta) 716 end do 717 718 dgxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r1 719 dgxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r3 720 d2gxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r2 721 d2gxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r4 722 d2gxdt(1,3,ilmn,ia,ispinor) = scale*buffer_r5 723 else 724 buffer_i1 = zero 725 buffer_i2 = zero 726 buffer_i3 = zero 727 buffer_i4 = zero 728 buffer_i5 = zero 729 do ipw=1,npw 730 aux_i = scali(ipw)*ffnl(ipw,1+igamma,ilmn) 731 buffer_i1 = buffer_i1 + aux_i 732 buffer_i2 = buffer_i2 + aux_i*kpg(ipw,ibeta) 733 aux_i = scali(ipw)*ffnl(ipw,1+idelta,ilmn) 734 buffer_i3 = buffer_i3 + aux_i 735 buffer_i4 = buffer_i4 + aux_i*kpg(ipw,ibeta) 736 aux_i = scali(ipw)*ffnl(ipw,4+idelgam,ilmn) 737 buffer_i5 = buffer_i5 + aux_i*kpg(ipw,ibeta) 738 end do 739 740 dgxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i1 741 dgxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i3 742 d2gxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i2 743 d2gxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i4 744 d2gxdt(1,3,ilmn,ia,ispinor) =-scale*buffer_i5 745 end if 746 end if 747 end do 748 end if 749 750 ! -------------------------------------------------------------------- 751 ! CHOICE= 4, 24 -- SIGNS= 1 752 ! Accumulate d2Gxdt --- 2nd derivative wrt 2 atm pos. --- for all dirs 753 ! -------------------------------------------------------------------- 754 if ((signs==1).and.(choice_==4.or.choice_==24)) then 755 mu0=1 756 do ilmn=1,nlmn 757 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 758 scale=wt;if (il>1) scale=-scale 759 scale2=two_pi2*scale 760 if (cplex==2) then 761 buffer_ra = zero ; buffer_rb = zero 762 buffer_rc = zero ; buffer_rd = zero 763 buffer_re = zero ; buffer_rf = zero 764 buffer_ia = zero ; buffer_ib = zero 765 buffer_ic = zero ; buffer_id = zero 766 buffer_ie = zero ; buffer_if = zero 767 do ipw=1,npw 768 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 769 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 770 buffer_ra = buffer_ra - aux_r*kpg(ipw,4) 771 buffer_rb = buffer_rb - aux_r*kpg(ipw,5) 772 buffer_rc = buffer_rc - aux_r*kpg(ipw,6) 773 buffer_rd = buffer_rd - aux_r*kpg(ipw,7) 774 buffer_re = buffer_re - aux_r*kpg(ipw,8) 775 buffer_rf = buffer_rf - aux_r*kpg(ipw,9) 776 buffer_ia = buffer_ia - aux_i*kpg(ipw,4) 777 buffer_ib = buffer_ib - aux_i*kpg(ipw,5) 778 buffer_ic = buffer_ic - aux_i*kpg(ipw,6) 779 buffer_id = buffer_id - aux_i*kpg(ipw,7) 780 buffer_ie = buffer_ie - aux_i*kpg(ipw,8) 781 buffer_if = buffer_if - aux_i*kpg(ipw,9) 782 end do 783 if (parity) then 784 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ra 785 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ia 786 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb 787 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_ib 788 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc 789 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_ic 790 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd 791 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_id 792 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re 793 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_ie 794 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf 795 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_if 796 else 797 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_ia 798 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ra 799 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib 800 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb 801 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic 802 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc 803 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id 804 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd 805 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie 806 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re 807 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if 808 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf 809 end if 810 else 811 if (parity) then 812 buffer_ra = zero ; buffer_rb = zero 813 buffer_rc = zero ; buffer_rd = zero 814 buffer_re = zero ; buffer_rf = zero 815 do ipw=1,npw 816 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 817 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 818 buffer_ra = buffer_ra - aux_r*kpg(ipw,4) 819 buffer_rb = buffer_rb - aux_r*kpg(ipw,5) 820 buffer_rc = buffer_rc - aux_r*kpg(ipw,6) 821 buffer_rd = buffer_rd - aux_r*kpg(ipw,7) 822 buffer_re = buffer_re - aux_r*kpg(ipw,8) 823 buffer_rf = buffer_rf - aux_r*kpg(ipw,9) 824 end do 825 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ra 826 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb 827 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc 828 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd 829 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re 830 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf 831 else 832 buffer_ia = zero ; buffer_ib = zero 833 buffer_ic = zero ; buffer_id = zero 834 buffer_ie = zero ; buffer_if = zero 835 do ipw=1,npw 836 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 837 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 838 buffer_ia = buffer_ia - aux_i*kpg(ipw,4) 839 buffer_ib = buffer_ib - aux_i*kpg(ipw,5) 840 buffer_ic = buffer_ic - aux_i*kpg(ipw,6) 841 buffer_id = buffer_id - aux_i*kpg(ipw,7) 842 buffer_ie = buffer_ie - aux_i*kpg(ipw,8) 843 buffer_if = buffer_if - aux_i*kpg(ipw,9) 844 end do 845 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_ia 846 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib 847 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic 848 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id 849 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie 850 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if 851 end if 852 end if 853 end do 854 end if 855 856 ! -------------------------------------------------------------------- 857 ! CHOICE= 5, 51, 52, 53, 54, 55, 8, 81 -- SIGNS= 1 858 ! Accumulate dGxdt --- derivative wrt k --- for all directions 859 ! -------------------------------------------------------------------- 860 if ((signs==1).and.& 861 & (choice_==5 .or.choice_==51.or.choice_==52.or.choice_==53.or.& 862 & choice_==54.or.choice_==55.or.choice_==8 .or.choice_==81)) then 863 mu0=1 864 if (choice_==54) mu0=4 865 if (choice_==55) mu0=7 866 do ilmn=1,nlmn 867 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 868 scale=wt;if (il>1) scale=-scale 869 if (cplex==2) then 870 buffer_r1 = zero ; buffer_r2 = zero 871 buffer_r3 = zero 872 buffer_i1 = zero ; buffer_i2 = zero 873 buffer_i3 = zero 874 do ipw=1,npw 875 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn) 876 buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn) 877 buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn) 878 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn) 879 buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn) 880 buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn) 881 end do 882 if (parity) then 883 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1 884 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 885 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2 886 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 887 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3 888 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 889 else 890 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_i1 891 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1 892 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2 893 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2 894 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3 895 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3 896 end if 897 else 898 cplex_dgxdt(mu0:mu0+2) = 2 ! Warning dgxdt is here pure imaginary 899 if (parity) then 900 buffer_i1 = zero ; buffer_i2 = zero 901 buffer_i3 = zero 902 do ipw=1,npw 903 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn) 904 buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn) 905 buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn) 906 end do 907 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 908 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 909 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 910 else 911 buffer_r1 = zero ; buffer_r2 = zero 912 buffer_r3 = zero 913 do ipw=1,npw 914 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn) 915 buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn) 916 buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn) 917 end do 918 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1 919 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2 920 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3 921 end if 922 end if 923 end do 924 end if 925 926 ! -------------------------------------------------------------------- 927 ! CHOICE= 5, 51, 52, 53, 54, 8, 81 -- SIGNS= 2 928 ! Accumulate dGxdt --- derivative wrt k --- 929 ! for direction IDIR (CHOICE=5,51,52) 930 ! for direction IDIR2 (CHOICE=54,81), only if choice > 0 931 ! for directions IDIR+1 & IDIR+2 mod 3 (CHOICE=53) 932 ! for 2 directions (CHOICE=8), only if choice > 0 933 ! -------------------------------------------------------------------- 934 if ((signs==2).and.& 935 & (choice_==5.or.choice_==51.or.choice_==52.or.choice_==53.or.choice_==54.or. & 936 & choice==8.or.choice==81)) then ! Note the use of choice and not choice_=abs(choice) 937 mu0=1 938 if (choice_==5.or.choice_==51.or.choice_==52) then 939 ! We compute the derivative in one direction 940 ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir 941 ndir=1 942 end if 943 if (choice_==53) then 944 ! Case 53 though need multiple directions and here ffnl will contain the 945 ! derivative information in locations 2,3, and 4 corresponding to idir = 1, 2, 946 ! and 3. Moreover, choice 53 needs the derivatives in direction idir+1 and idir+2 mod 3. 947 ! The parameter vector ffnl_dir_dat contains the necessary translations in 948 ! locations 1,2 for idir=1; 3,4 for idir=2; and 5,6 for idir=3. 949 ffnl_dir(1)=ffnl_dir_dat(2*idir-1);ffnl_dir(2)=ffnl_dir_dat(2*idir) 950 ndir=2 951 end if 952 if (choice_==54) then 953 ! We compute the derivative in one direction (the second one) 954 ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir2(idir) 955 ndir=1 956 end if 957 if (choice_==8) then 958 ! We compute the derivative in two directions 959 ! depending on the input idir argument 960 ffnl_dir(1)=idir1(idir);ffnl_dir(2)=idir2(idir) 961 ndir=2 962 end if 963 if (choice_==81) then 964 ! We compute the derivative in one direction (the second one) 965 ffnl_dir(1)=idir2(idir) 966 ndir=1 967 end if 968 do mua=1,ndir 969 mu=mu0+mua-1 970 do ilmn=1,nlmn 971 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 972 scale=wt;if (il>1) scale=-scale 973 if (cplex==2) then 974 buffer_r1 = zero ; buffer_i1 = zero 975 do ipw=1,npw 976 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 977 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 978 end do 979 if (parity) then 980 dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1 981 dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_i1 982 else 983 dgxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1 984 dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_r1 985 end if 986 else 987 cplex_dgxdt(mu) = 2 ! Warning dgxdt is here pure imaginary 988 if (parity) then 989 buffer_i1 = zero 990 do ipw=1,npw 991 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 992 end do 993 dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1 994 else 995 buffer_r1 = zero 996 do ipw=1,npw 997 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 998 end do 999 dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1 1000 end if 1001 end if 1002 end do 1003 end do 1004 end if 1005 1006 ! -------------------------------------------------------------------- 1007 ! CHOICE= 54 -- SIGNS= 1 1008 ! Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- for all dirs 1009 ! -------------------------------------------------------------------- 1010 if ((signs==1).and.(choice_==54)) then 1011 mu0=1 1012 do ilmn=1,nlmn 1013 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1014 scale=wt;if (il>1) scale=-scale 1015 scale2=two_pi*scale 1016 if (cplex==2) then 1017 do mua=1,3 ! atm. pos 1018 do mub=1,3 ! k 1019 mu=mu0+(mub-1)+3*(mua-1) 1020 buffer_r1 = zero ; buffer_i1 = zero 1021 do ipw=1,npw 1022 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1023 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1024 end do 1025 if (parity) then 1026 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 1027 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1028 else 1029 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_r1 1030 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 1031 end if 1032 end do 1033 end do 1034 else 1035 cplex_d2gxdt(mu0:mu0+8) = 2 ! Warning d2gxdt is here pure imaginary 1036 if (parity) then 1037 do mua=1,3 ! atm. pos 1038 do mub=1,3 ! k 1039 mu=mu0+(mub-1)+3*(mua-1) 1040 buffer_r1 = zero 1041 do ipw=1,npw 1042 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1043 end do 1044 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1045 end do 1046 end do 1047 else 1048 do mua=1,3 ! atm. pos 1049 do mub=1,3 ! k 1050 mu=mu0+(mub-1)+3*(mua-1) 1051 buffer_i1 = zero 1052 do ipw=1,npw 1053 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1054 end do 1055 d2gxdt(1,mu,ilmn,ia,ispinor) = -scale2*buffer_i1 1056 end do 1057 end do 1058 end if 1059 end if 1060 end do 1061 end if 1062 1063 ! -------------------------------------------------------------------- 1064 ! CHOICE= 54 -- SIGNS= 2 1065 ! Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- for direction IDIR 1066 ! -------------------------------------------------------------------- 1067 if ((signs==2).and.(choice_==54)) then 1068 mu0=1 1069 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1070 mua=idir1(idir) ! atm. pos 1071 mub=1;if(dimffnl>2) mub=idir2(idir) ! k 1072 do ilmn=1,nlmn 1073 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1074 scale=wt;if (il>1) scale=-scale 1075 scale2=two_pi*scale 1076 if (cplex==2) then 1077 buffer_r1 = zero ; buffer_i1 = zero 1078 do ipw=1,npw 1079 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1080 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1081 end do 1082 if (parity) then 1083 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 1084 d2gxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 1085 else 1086 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_r1 1087 d2gxdt(2,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 1088 end if 1089 else 1090 cplex_d2gxdt(mu0) = 2 ! Warning d2gxdt is here pure imaginary 1091 if (parity) then 1092 buffer_r1 = zero 1093 do ipw=1,npw 1094 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1095 end do 1096 d2gxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 1097 else 1098 buffer_i1 = zero 1099 do ipw=1,npw 1100 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 1101 end do 1102 d2gxdt(1,mu0,ilmn,ia,ispinor) = -scale2*buffer_i1 1103 end if 1104 end if 1105 end do 1106 end if 1107 1108 ! -------------------------------------------------------------------- 1109 ! CHOICE= 55 -- SIGNS= 1 1110 ! Accumulate d2Gxdt --- 2nd derivative wrt k and strains --- for all dirs 1111 ! -------------------------------------------------------------------- 1112 if ((signs==1).and.(choice_==55)) then 1113 mu0=1 1114 do ilmn=1,nlmn 1115 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1116 scale=wt;if (il>1) scale=-scale 1117 if (cplex==2) then 1118 do mua=1,6 ! strain 1119 do mub=1,3 ! k 1120 mu=mu0+(mub-1)+3*(mua-1) 1121 buffer_r1 = zero ; buffer_i1 = zero 1122 do ipw=1,npw 1123 aux_r=ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub) 1124 buffer_r1 = buffer_r1 + aux_r*scalr(ipw) 1125 buffer_i1 = buffer_i1 + aux_r*scali(ipw) 1126 end do 1127 if (parity) then 1128 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1 1129 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_i1 1130 else 1131 d2gxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1 1132 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_r1 1133 end if 1134 end do 1135 end do 1136 else 1137 cplex_d2gxdt(mu0:mu0+17) = 2 ! Warning d2gxdt is here pure imaginary 1138 if (parity) then 1139 do mua=1,6 ! strain 1140 do mub=1,3 ! k 1141 mu=mu0+(mub-1)+3*(mua-1) 1142 buffer_i1 = zero 1143 do ipw=1,npw 1144 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+mua,ilmn)*kpg(ipw,mub) 1145 end do 1146 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1 1147 end do 1148 end do 1149 else 1150 do mua=1,6 ! strain 1151 do mub=1,3 ! k 1152 mu=mu0+(mub-1)+3*(mua-1) 1153 buffer_r1 = zero 1154 do ipw=1,npw 1155 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+mua,ilmn)*kpg(ipw,mub) 1156 end do 1157 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1 1158 end do 1159 end do 1160 end if 1161 end if 1162 end do 1163 end if 1164 1165 ! -------------------------------------------------------------------- 1166 ! CHOICE= 6 -- SIGNS= 1 1167 ! Accumulate d2Gxdt --- 2nd derivative wrt atm. pos and strain --- for all dirs 1168 ! -------------------------------------------------------------------- 1169 if ((signs==1).and.(choice_==6)) then 1170 mu0=1;if (choice_==6) mu0=37 1171 ABI_MALLOC(scalarr,(npw,4)) 1172 ABI_MALLOC(scalari,(npw,4)) 1173 do ilmn=1,nlmn 1174 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1175 scale=wt;if (il>1) scale=-scale 1176 scale2=pi*scale 1177 do mu=1,4 1178 do ipw=1,npw 1179 scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn) 1180 scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn) 1181 end do 1182 end do 1183 if(cplex==2) then 1184 do mub=1,6 1185 nub1=alpha(mub);nub2=beta(mub) 1186 do mua=1,3 1187 mu=mu0+(mua-1)+3*(mub-1) 1188 buffer_r1 = zero ; buffer_i1 = zero 1189 do ipw=1,npw 1190 buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + & 1191 & scalarr(ipw,1+nub2)*kpg(ipw,nub1)) 1192 buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + & 1193 & scalari(ipw,1+nub2)*kpg(ipw,nub1)) 1194 end do 1195 if (parity) then 1196 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1 1197 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_r1 1198 else 1199 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1200 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1 1201 end if 1202 end do 1203 end do 1204 else 1205 if (parity) then 1206 do mub=1,6 1207 nub1=alpha(mub);nub2=beta(mub) 1208 do mua=1,3 1209 mu=mu0+(mua-1)+3*(mub-1) 1210 buffer_i1 = zero 1211 do ipw=1,npw 1212 buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + & 1213 & scalari(ipw,1+nub2)*kpg(ipw,nub1)) 1214 end do 1215 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1 1216 end do 1217 end do 1218 else 1219 do mub=1,6 1220 nub1=alpha(mub);nub2=beta(mub) 1221 do mua=1,3 1222 mu=mu0+(mua-1)+3*(mub-1) 1223 buffer_r1 = zero 1224 do ipw=1,npw 1225 buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + & 1226 & scalarr(ipw,1+nub2)*kpg(ipw,nub1)) 1227 end do 1228 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1229 end do 1230 end do 1231 end if 1232 end if 1233 end do 1234 ABI_FREE(scalarr) 1235 ABI_FREE(scalari) 1236 end if 1237 1238 ! -------------------------------------------------------------------- 1239 ! CHOICE= 6 -- SIGNS= 1 1240 ! Accumulate d2Gxdt --- 2nd derivative wrt 2 strains --- for all dirs 1241 ! -------------------------------------------------------------------- 1242 if ((signs==1).and.(choice_==6)) then 1243 mu0=1 1244 ABI_MALLOC(scalarr,(npw,10)) 1245 ABI_MALLOC(scalari,(npw,10)) 1246 do ilmn=1,nlmn 1247 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1248 scale=wt;if (il>1) scale=-scale 1249 scale2=quarter*scale 1250 do mu=1,10 1251 do ipw=1,npw 1252 scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn) 1253 scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn) 1254 end do 1255 end do 1256 if(cplex==2) then 1257 do mub=1,6 1258 nub1=alpha(mub);nub2=beta(mub) 1259 do mua=1,6 1260 mu=mu0+(mua-1)+6*(mub-1) 1261 nua1=alpha(mua);nua2=beta(mua) 1262 gama=gamma(nub1,nua1) 1263 gamb=gamma(nub2,nua1) 1264 gamc=gamma(nub1,nua2) 1265 gamd=gamma(nub2,nua2) 1266 buffer_r1 = zero ; buffer_i1 = zero 1267 do ipw=1,npw 1268 kpga=kpg(ipw,nub1) 1269 kpgb=kpg(ipw,nub2) 1270 kpgc=kpg(ipw,nua1) 1271 kpgd=kpg(ipw,nua2) 1272 buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd & 1273 & + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc 1274 buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd & 1275 & + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc 1276 end do 1277 if (parity) then 1278 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1279 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1 1280 else 1281 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 1282 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1283 end if 1284 end do 1285 end do 1286 else 1287 if (parity) then 1288 do mub=1,6 1289 nub1=alpha(mub);nub2=beta(mub) 1290 do mua=1,6 1291 mu=mu0+(mua-1)+6*(mub-1) 1292 nua1=alpha(mua);nua2=beta(mua) 1293 gama=gamma(nub1,nua1) 1294 gamb=gamma(nub2,nua1) 1295 gamc=gamma(nub1,nua2) 1296 gamd=gamma(nub2,nua2) 1297 buffer_r1 = zero 1298 do ipw=1,npw 1299 kpga=kpg(ipw,nub1) 1300 kpgb=kpg(ipw,nub2) 1301 kpgc=kpg(ipw,nua1) 1302 kpgd=kpg(ipw,nua2) 1303 buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd & 1304 & + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc 1305 end do 1306 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 1307 end do 1308 end do 1309 else 1310 do mub=1,6 1311 nub1=alpha(mub);nub2=beta(mub) 1312 do mua=1,6 1313 mu=mu0+(mua-1)+6*(mub-1) 1314 nua1=alpha(mua);nua2=beta(mua) 1315 gama=gamma(nub1,nua1) 1316 gamb=gamma(nub2,nua1) 1317 gamc=gamma(nub1,nua2) 1318 gamd=gamma(nub2,nua2) 1319 buffer_i1 = zero 1320 do ipw=1,npw 1321 kpga=kpg(ipw,nub1) 1322 kpgb=kpg(ipw,nub2) 1323 kpgc=kpg(ipw,nua1) 1324 kpgd=kpg(ipw,nua2) 1325 buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd & 1326 & + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc 1327 end do 1328 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 1329 end do 1330 end do 1331 end if 1332 end if 1333 end do 1334 ABI_FREE(scalarr) 1335 ABI_FREE(scalari) 1336 end if 1337 1338 ! -------------------------------------------------------------------- 1339 ! CHOICE= 8, 81 -- SIGNS= 1 1340 ! Accumulate d2Gxdt --- 2nd derivative wrt 2 k wave vect. --- for all dirs 1341 ! -------------------------------------------------------------------- 1342 if ((signs==1).and.(choice_==8.or.choice_==81)) then 1343 mu0=1 1344 do ilmn=1,nlmn 1345 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1346 scale=wt;if (il>1) scale=-scale 1347 if (cplex==2) then 1348 buffer_ra = zero ; buffer_rb = zero 1349 buffer_rc = zero ; buffer_rd = zero 1350 buffer_re = zero ; buffer_rf = zero 1351 buffer_ia = zero ; buffer_ib = zero 1352 buffer_ic = zero ; buffer_id = zero 1353 buffer_ie = zero ; buffer_if = zero 1354 do ipw=1,npw 1355 buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn) 1356 buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn) 1357 buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn) 1358 buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn) 1359 buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn) 1360 buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn) 1361 buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn) 1362 buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn) 1363 buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn) 1364 buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn) 1365 buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn) 1366 buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn) 1367 end do 1368 if (parity) then 1369 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_ra 1370 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_ia 1371 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb 1372 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_ib 1373 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc 1374 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_ic 1375 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd 1376 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_id 1377 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re 1378 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_ie 1379 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf 1380 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_if 1381 else 1382 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_ia 1383 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_ra 1384 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib 1385 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb 1386 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic 1387 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc 1388 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id 1389 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd 1390 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie 1391 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_re 1392 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if 1393 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf 1394 end if 1395 else 1396 if (parity) then 1397 buffer_ra = zero ; buffer_rb = zero 1398 buffer_rc = zero ; buffer_rd = zero 1399 buffer_re = zero ; buffer_rf = zero 1400 do ipw=1,npw 1401 buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn) 1402 buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn) 1403 buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn) 1404 buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn) 1405 buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn) 1406 buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn) 1407 end do 1408 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_ra 1409 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb 1410 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc 1411 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd 1412 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re 1413 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf 1414 else 1415 buffer_ia = zero ; buffer_ib = zero 1416 buffer_ic = zero ; buffer_id = zero 1417 buffer_ie = zero ; buffer_if = zero 1418 do ipw=1,npw 1419 buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn) 1420 buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn) 1421 buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn) 1422 buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn) 1423 buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn) 1424 buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn) 1425 end do 1426 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_ia 1427 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib 1428 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic 1429 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id 1430 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie 1431 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if 1432 end if 1433 end if 1434 end do 1435 end if 1436 1437 ! -------------------------------------------------------------------- 1438 ! CHOICE= 8, 81 -- SIGNS= 2 1439 ! Accumulate d2Gxdt --- 2nd derivative wrt k --- for direction IDIR (Voigt not.) 1440 ! -------------------------------------------------------------------- 1441 if ((signs==2).and.(choice_==8.or.choice_==81)) then 1442 mu0=1 1443 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 1444 mua=idir1(idir);mub=idir2(idir) 1445 !idir->(Voigt) (11->1,22->2,33->3,32->4,31->5,21->6) 1446 ffnl_dir(1)=gamma(mua,mub) 1447 do ilmn=1,nlmn 1448 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1449 scale=wt;if (il>1) scale=-scale 1450 if (cplex==2) then 1451 buffer_r1 = zero ; buffer_i1 = zero 1452 do ipw=1,npw 1453 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 1454 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 1455 end do 1456 if (parity) then 1457 d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1458 d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1 1459 else 1460 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 1461 d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1462 end if 1463 else 1464 if (parity) then 1465 buffer_r1 = zero 1466 do ipw=1,npw 1467 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 1468 end do 1469 d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1470 else 1471 buffer_i1 = zero 1472 do ipw=1,npw 1473 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 1474 end do 1475 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 1476 end if 1477 end if 1478 end do 1479 end if 1480 1481 ! -------------------------------------------------------------------- 1482 ! END CHOICES 1483 ! -------------------------------------------------------------------- 1484 1485 end do ! End loop on atoms 1486 1487 ! Deallocate temporary space 1488 ABI_FREE(scali) 1489 ABI_FREE(scalr) 1490 1491 end do ! End loop on spinorial components 1492 1493 1494 ! ========================================================================== 1495 ! ========== OPENMP VERSION ================================================ 1496 ! ========================================================================== 1497 ! Parallelization OpenMP on npw loops 1498 else 1499 1500 !$OMP PARALLEL PRIVATE(il,ilmn,ipw,jpw,parity,ndir,ffnl_dir), & 1501 !$OMP PRIVATE(mu,mua,mub,nua1,nua2,nub1,nub2), & 1502 !$OMP PRIVATE(kpga,kpgb,kpgc,kpgd), & 1503 !$OMP PRIVATE(scale,scale2,mu0), & 1504 !$OMP PRIVATE(ispinor,ipwshft,ia,iaph3d) 1505 1506 ! Loop on spinorial components 1507 do ispinor =1,nspinor 1508 ipwshft=(ispinor-1)*npw 1509 1510 ! Allocate work space 1511 !$OMP SECTIONS 1512 !$OMP SECTION 1513 ABI_MALLOC(scali,(npw)) 1514 !$OMP SECTION 1515 ABI_MALLOC(scalr,(npw)) 1516 !$OMP END SECTIONS 1517 1518 ! Loop on atoms (blocking) 1519 do ia=1,nincat 1520 iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1 1521 1522 ! Compute Sum_g[c(g).exp(2pi.i.g.R)] 1523 !$OMP DO 1524 do ipw=ipw0,npw 1525 jpw=ipw+ipwshft 1526 scalr(ipw)=(vect(1,jpw)*ph3d(1,ipw,iaph3d)-vect(2,jpw)*ph3d(2,ipw,iaph3d)) 1527 scali(ipw)=(vect(2,jpw)*ph3d(1,ipw,iaph3d)+vect(1,jpw)*ph3d(2,ipw,iaph3d)) 1528 end do 1529 !$OMP END DO 1530 !$OMP SINGLE 1531 if (ipw0==2) then 1532 scalr(1)=half*vect(1,1+ipwshft)*ph3d(1,1,iaph3d) 1533 scali(1)=half*vect(1,1+ipwshft)*ph3d(2,1,iaph3d) 1534 end if 1535 !$OMP END SINGLE 1536 1537 ! -------------------------------------------------------------------- 1538 ! ALL CHOICES 1539 ! Accumulate Gx 1540 ! -------------------------------------------------------------------- 1541 if (choice>=0) then ! JWZ to check: I dont think 53 needs this 1542 do ilmn=1,nlmn 1543 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1544 scale=wt;if (il>1) scale=-scale 1545 if (cplex==2) then 1546 !$OMP SINGLE 1547 buffer_r1 = zero ; buffer_i1 = zero 1548 !$OMP END SINGLE 1549 !$OMP DO & 1550 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 1551 do ipw=1,npw 1552 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn) 1553 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn) 1554 end do 1555 !$OMP END DO 1556 !$OMP SINGLE 1557 if (parity) then 1558 gx(1,ilmn,ia,ispinor) = scale*buffer_r1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_i1 1559 else 1560 gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_r1 1561 end if 1562 !$OMP END SINGLE 1563 else 1564 if (parity) then 1565 !$OMP SINGLE 1566 buffer_r1 = zero 1567 !$OMP END SINGLE 1568 !$OMP DO & 1569 !$OMP REDUCTION(+:buffer_r1) 1570 do ipw=1,npw 1571 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn) 1572 end do 1573 !$OMP END DO 1574 !$OMP SINGLE 1575 gx(1,ilmn,ia,ispinor) = scale*buffer_r1 1576 !$OMP END SINGLE 1577 else 1578 !$OMP SINGLE 1579 buffer_i1 = zero 1580 !$OMP END SINGLE 1581 !$OMP DO & 1582 !$OMP REDUCTION(+:buffer_i1) 1583 do ipw=1,npw 1584 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn) 1585 end do 1586 !$OMP END DO 1587 !$OMP SINGLE 1588 gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 1589 !$OMP END SINGLE 1590 end if 1591 end if 1592 end do 1593 end if 1594 1595 ! -------------------------------------------------------------------- 1596 ! CHOICE= 2, 4, 6, 23, 24, 54 -- SIGNS= 1 1597 ! Accumulate dGxdt --- derivative wrt atm pos. --- for all directions 1598 ! -------------------------------------------------------------------- 1599 if ((signs==1).and. & 1600 & (choice_==2.or.choice_==23.or.choice_==24.or.choice_==4.or.choice_==54.or.choice_==6)) then 1601 mu0 =1;if (choice_==23.or.choice_==6) mu0=7 1602 do ilmn=1,nlmn 1603 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1604 scale=wt;if (il>1) scale=-scale 1605 scale2=two_pi*scale 1606 if (cplex==2) then 1607 !$OMP SINGLE 1608 buffer_r1 = zero ; buffer_r2 = zero 1609 buffer_r3 = zero ; buffer_i1 = zero 1610 buffer_i2 = zero ; buffer_i3 = zero 1611 !$OMP END SINGLE 1612 !$OMP DO & 1613 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_i1,buffer_i2,buffer_i3),& 1614 !$OMP PRIVATE(aux_r,aux_i) 1615 do ipw=1,npw 1616 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 1617 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 1618 buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1) 1619 buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2) 1620 buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3) 1621 buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1) 1622 buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2) 1623 buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3) 1624 end do 1625 !$OMP END DO 1626 !$OMP SINGLE 1627 if (parity) then 1628 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1 1629 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_i1 1630 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2 1631 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_i2 1632 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3 1633 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_i3 1634 else 1635 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_i1 1636 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1 1637 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2 1638 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2 1639 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3 1640 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3 1641 end if 1642 !$OMP END SINGLE 1643 else 1644 if (parity) then 1645 !$OMP SINGLE 1646 buffer_r1 = zero ; buffer_r2 = zero ; buffer_r3 = zero 1647 !$OMP END SINGLE 1648 !$OMP DO & 1649 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3), & 1650 !$OMP PRIVATE(aux_r,aux_i) 1651 do ipw=1,npw 1652 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 1653 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 1654 buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1) 1655 buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2) 1656 buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3) 1657 end do 1658 !$OMP END DO 1659 !$OMP SINGLE 1660 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1 1661 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2 1662 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3 1663 !$OMP END SINGLE 1664 else 1665 !$OMP SINGLE 1666 buffer_i1 = zero ; buffer_i2 = zero ; buffer_i3 = zero 1667 !$OMP END SINGLE 1668 !$OMP DO & 1669 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3), & 1670 !$OMP PRIVATE(aux_r,aux_i) 1671 do ipw=1,npw 1672 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 1673 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 1674 buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1) 1675 buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2) 1676 buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3) 1677 end do 1678 !$OMP SINGLE 1679 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_i1 1680 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2 1681 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3 1682 !$OMP END SINGLE 1683 end if 1684 end if 1685 end do 1686 end if 1687 1688 ! -------------------------------------------------------------------- 1689 ! CHOICE= 2 -- SIGNS= 2 1690 ! Accumulate dGxdt --- derivative wrt atm pos. --- for direction IDIR 1691 ! -------------------------------------------------------------------- 1692 if ((signs==2).and.(choice_==2)) then 1693 mu0=1 1694 do ilmn=1,nlmn 1695 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1696 scale=wt;if (il>1) scale=-scale 1697 scale2=two_pi*scale 1698 if (cplex==2) then 1699 !$OMP SINGLE 1700 buffer_r1 = zero ; buffer_i1 = zero 1701 !$OMP END SINGLE 1702 !$OMP DO & 1703 !$OMP REDUCTION(+:buffer_r1,buffer_i1), & 1704 !$OMP PRIVATE(aux_r,aux_i) 1705 do ipw=1,npw 1706 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 1707 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 1708 buffer_r1 = buffer_r1 - aux_i*kpg(ipw,idir) 1709 buffer_i1 = buffer_i1 + aux_r*kpg(ipw,idir) 1710 end do 1711 !$OMP END DO 1712 !$OMP SINGLE 1713 if (parity) then 1714 dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 1715 dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_i1 1716 else 1717 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 1718 dgxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 1719 end if 1720 !$OMP END SINGLE 1721 else 1722 if (parity) then 1723 !$OMP SINGLE 1724 buffer_r1 = zero 1725 !$OMP END SINGLE 1726 !$OMP DO & 1727 !$OMP REDUCTION(+:buffer_r1) 1728 do ipw=1,npw 1729 buffer_r1 = buffer_r1 - scali(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir) 1730 end do 1731 !$OMP END DO 1732 !$OMP SINGLE 1733 dgxdt(1,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 1734 !$OMP END SINGLE 1735 else 1736 !$OMP SINGLE 1737 buffer_i1 = zero 1738 !$OMP END SINGLE 1739 !$OMP DO & 1740 !$OMP REDUCTION(+:buffer_i1) 1741 do ipw=1,npw 1742 buffer_i1 = buffer_i1 + scalr(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,idir) 1743 end do 1744 !$OMP END DO 1745 !$OMP SINGLE 1746 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 1747 !$OMP END SINGLE 1748 end if 1749 end if 1750 end do 1751 end if 1752 1753 ! -------------------------------------------------------------------- 1754 ! CHOICE= 22 and 25 -- SIGNS= 2 1755 ! Accumulate dGxdt --- derivative wrt atm pos. and one or two q-vectors--- 1756 ! for directions fixed by IDIR and qdir 1757 ! -------------------------------------------------------------------- 1758 if (signs==2.and.(choice_==22.or.choice_==25)) then 1759 !Use same notation as the notes for clarity 1760 if(choice_==22) ialpha=idir 1761 if(choice_==25) ialpha=nalpha(idir) 1762 mu0=1 1763 do ilmn=1,nlmn 1764 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1765 scale=wt;if (il>1) scale=-scale 1766 if (cplex==2) then 1767 !$OMP SINGLE 1768 buffer_r1 = zero ; buffer_i1 = zero 1769 !$OMP END SINGLE 1770 !$OMP DO & 1771 !$OMP REDUCTION(+:buffer_r1,buffer_i1), & 1772 !$OMP PRIVATE(aux_r,aux_i) 1773 do ipw=1,npw 1774 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 1775 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 1776 buffer_r1 = buffer_r1 + aux_r*kpg(ipw,ialpha) 1777 buffer_i1 = buffer_i1 + aux_i*kpg(ipw,ialpha) 1778 end do 1779 !$OMP END DO 1780 !$OMP SINGLE 1781 if (parity) then 1782 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1783 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1 1784 else 1785 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 1786 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1787 end if 1788 !$OMP END SINGLE 1789 else 1790 if (parity) then 1791 !$OMP SINGLE 1792 buffer_r1 = zero 1793 !$OMP END SINGLE 1794 !$OMP DO & 1795 !$OMP REDUCTION(+:buffer_r1) 1796 do ipw=1,npw 1797 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,ialpha) 1798 end do 1799 !$OMP END DO 1800 !$OMP SINGLE 1801 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1802 !$OMP END SINGLE 1803 else 1804 !$OMP SINGLE 1805 buffer_i1 = zero 1806 !$OMP END SINGLE 1807 !$OMP DO & 1808 !$OMP REDUCTION(+:buffer_i1) 1809 do ipw=1,npw 1810 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)*kpg(ipw,ialpha) 1811 end do 1812 !$OMP END DO 1813 !$OMP SINGLE 1814 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 1815 !$OMP END SINGLE 1816 end if 1817 end if 1818 end do 1819 end if 1820 1821 ! -------------------------------------------------------------------- 1822 ! CHOICE= 3, 6, 23, 55 -- SIGNS= 1 1823 ! Accumulate dGxdt --- derivative wrt strain --- for all directions 1824 ! -------------------------------------------------------------------- 1825 if ((signs==1).and.(choice_==3.or.choice_==6.or.choice_==23.or.choice_==55)) then 1826 mu0=1 1827 do ilmn=1,nlmn 1828 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1829 scale=wt;if (il>1) scale=-scale 1830 scale2=half*scale 1831 if (cplex==2) then 1832 !$OMP SINGLE 1833 buffer_r1 = zero ; buffer_r2 = zero 1834 buffer_r3 = zero ; buffer_r4 = zero 1835 buffer_r5 = zero ; buffer_r6 = zero 1836 buffer_i1 = zero ; buffer_i2 = zero 1837 buffer_i3 = zero ; buffer_i4 = zero 1838 buffer_i5 = zero ; buffer_i6 = zero 1839 !$OMP END SINGLE 1840 !$OMP DO & 1841 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5,buffer_r6), & 1842 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5,buffer_i6), & 1843 !$OMP PRIVATE(aux_r2,aux_r3,aux_r4,aux_i2,aux_i3,aux_i4) 1844 do ipw=1,npw 1845 aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn) 1846 aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn) 1847 aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn) 1848 aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn) 1849 aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn) 1850 aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn) 1851 buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1) 1852 buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2) 1853 buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3) 1854 buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3) 1855 buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3) 1856 buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2) 1857 buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1) 1858 buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2) 1859 buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3) 1860 buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3) 1861 buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3) 1862 buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2) 1863 end do 1864 !$OMP END DO 1865 !$OMP SINGLE 1866 if (parity) then 1867 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1 1868 dgxdt(2,mu0 ,ilmn,ia,ispinor) =-scale*buffer_i1 1869 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2 1870 dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2 1871 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3 1872 dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3 1873 dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4 1874 dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_i4 1875 dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5 1876 dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_i5 1877 dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6 1878 dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_i6 1879 else 1880 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 1881 dgxdt(2,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1 1882 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 1883 dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2 1884 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 1885 dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3 1886 dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4 1887 dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4 1888 dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5 1889 dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5 1890 dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6 1891 dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6 1892 end if 1893 !$OMP END SINGLE 1894 else 1895 if (parity) then 1896 !$OMP SINGLE 1897 buffer_r1 = zero ; buffer_r2 = zero 1898 buffer_r3 = zero ; buffer_r4 = zero 1899 buffer_r5 = zero ; buffer_r6 = zero 1900 !$OMP END SINGLE 1901 !$OMP DO & 1902 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5,buffer_r6), & 1903 !$OMP PRIVATE(aux_r2,aux_r3,aux_r4) 1904 do ipw=1,npw 1905 aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn) 1906 aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn) 1907 aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn) 1908 buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1) 1909 buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2) 1910 buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3) 1911 buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3) 1912 buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3) 1913 buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2) 1914 end do 1915 !$OMP SINGLE 1916 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1 1917 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2 1918 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3 1919 dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4 1920 dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5 1921 dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6 1922 !$OMP END SINGLE 1923 else 1924 !$OMP SINGLE 1925 buffer_i1 = zero ; buffer_i2 = zero 1926 buffer_i3 = zero ; buffer_i4 = zero 1927 buffer_i5 = zero ; buffer_i6 = zero 1928 !$OMP END SINGLE 1929 !$OMP DO & 1930 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5,buffer_i6), & 1931 !$OMP PRIVATE(aux_i2,aux_i3,aux_i4) 1932 do ipw=1,npw 1933 aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn) 1934 aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn) 1935 aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn) 1936 buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1) 1937 buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2) 1938 buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3) 1939 buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3) 1940 buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3) 1941 buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2) 1942 end do 1943 !$OMP END DO 1944 !$OMP SINGLE 1945 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 1946 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 1947 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 1948 dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4 1949 dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5 1950 dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6 1951 !$OMP END SINGLE 1952 end if 1953 end if 1954 end do 1955 end if 1956 1957 ! -------------------------------------------------------------------- 1958 ! CHOICE= 3 -- SIGNS= 2 1959 ! Accumulate dGxdt --- derivative wrt strain --- for direction IDIR 1960 ! -------------------------------------------------------------------- 1961 if ((signs==2).and.(choice_==3)) then 1962 mu0=1;ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir 1963 do ilmn=1,nlmn 1964 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 1965 scale=wt;if (il>1) scale=-scale 1966 if (cplex==2) then 1967 !$OMP SINGLE 1968 buffer_r1 = zero ; buffer_i1 = zero 1969 !$OMP END SINGLE 1970 !$OMP DO & 1971 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 1972 do ipw=1,npw 1973 buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 1974 buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 1975 end do 1976 !$OMP END DO 1977 !$OMP SINGLE 1978 if (parity) then 1979 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1980 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1 1981 else 1982 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 1983 dgxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1984 end if 1985 !$OMP END SINGLE 1986 else 1987 if (parity) then 1988 !$OMP SINGLE 1989 buffer_r1 = zero 1990 !$OMP END SINGLE 1991 !$OMP DO & 1992 !$OMP REDUCTION(+:buffer_r1) 1993 do ipw=1,npw 1994 buffer_r1 = buffer_r1 - scalr(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 1995 end do 1996 !$OMP END DO 1997 !$OMP SINGLE 1998 dgxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 1999 !$OMP END SINGLE 2000 else 2001 !$OMP SINGLE 2002 buffer_i1 = zero 2003 !$OMP END SINGLE 2004 !$OMP DO & 2005 !$OMP REDUCTION(+:buffer_i1) 2006 do ipw=1,npw 2007 buffer_i1 = buffer_i1 - scali(ipw)*ffnl(ipw,1+ffnl_dir(1),ilmn) 2008 end do 2009 !$OMP END DO 2010 !$OMP SINGLE 2011 dgxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 2012 !$OMP END SINGLE 2013 end if 2014 end if 2015 end do 2016 end if 2017 2018 ! -------------------------------------------------------------------- 2019 ! CHOICE= 33 -- SIGNS= 2 2020 ! Accumulate dGxdt and d2Gxdt --- mixed derivative wrt strain and q-vector --- 2021 ! for directions fixed by idir and qdir 2022 ! -------------------------------------------------------------------- 2023 if ((signs==2).and.(choice_==33)) then 2024 !Use same notation as the notes for clarity 2025 ibeta=nalpha(idir) 2026 idelta=nbeta(idir) 2027 igamma=qdir 2028 idelgam=gamma(idelta,igamma) 2029 do ilmn=1,nlmn 2030 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2031 scale=wt;if (il>1) scale=-scale 2032 if (cplex==2) then 2033 !$OMP SINGLE 2034 buffer_r1 = zero ; buffer_i1 = zero 2035 buffer_r2 = zero ; buffer_i2 = zero 2036 buffer_r3 = zero ; buffer_i3 = zero 2037 buffer_r4 = zero ; buffer_i4 = zero 2038 buffer_r5 = zero ; buffer_i5 = zero 2039 !$OMP END SINGLE 2040 !$OMP DO & 2041 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5), & 2042 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5), & 2043 !$OMP PRIVATE(aux_r,aux_i) 2044 do ipw=1,npw 2045 aux_r = scalr(ipw)*ffnl(ipw,1+igamma,ilmn) 2046 aux_i = scali(ipw)*ffnl(ipw,1+igamma,ilmn) 2047 buffer_r1 = buffer_r1 + aux_r 2048 buffer_i1 = buffer_i1 + aux_i 2049 buffer_r2 = buffer_r2 + aux_r*kpg(ipw,ibeta) 2050 buffer_i2 = buffer_i2 + aux_i*kpg(ipw,ibeta) 2051 aux_r = scalr(ipw)*ffnl(ipw,1+idelta,ilmn) 2052 aux_i = scali(ipw)*ffnl(ipw,1+idelta,ilmn) 2053 buffer_r3 = buffer_r3 + aux_r 2054 buffer_i3 = buffer_i3 + aux_i 2055 buffer_r4 = buffer_r4 + aux_r*kpg(ipw,ibeta) 2056 buffer_i4 = buffer_i4 + aux_i*kpg(ipw,ibeta) 2057 aux_r = scalr(ipw)*ffnl(ipw,4+idelgam,ilmn) 2058 aux_i = scali(ipw)*ffnl(ipw,4+idelgam,ilmn) 2059 buffer_r5 = buffer_r5 + aux_r*kpg(ipw,ibeta) 2060 buffer_i5 = buffer_i5 + aux_i*kpg(ipw,ibeta) 2061 end do 2062 !$OMP END DO 2063 !$OMP SINGLE 2064 if (parity) then 2065 dgxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r1 2066 dgxdt(2,1,ilmn,ia,ispinor) = scale*buffer_i1 2067 dgxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r3 2068 dgxdt(2,2,ilmn,ia,ispinor) = scale*buffer_i3 2069 d2gxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r2 2070 d2gxdt(2,1,ilmn,ia,ispinor) = scale*buffer_i2 2071 d2gxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r4 2072 d2gxdt(2,2,ilmn,ia,ispinor) = scale*buffer_i4 2073 d2gxdt(1,3,ilmn,ia,ispinor) = scale*buffer_r5 2074 d2gxdt(2,3,ilmn,ia,ispinor) = scale*buffer_i5 2075 else 2076 dgxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i1 2077 dgxdt(2,1,ilmn,ia,ispinor) = scale*buffer_r1 2078 dgxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i3 2079 dgxdt(2,2,ilmn,ia,ispinor) = scale*buffer_r3 2080 d2gxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i2 2081 d2gxdt(2,1,ilmn,ia,ispinor) = scale*buffer_r2 2082 d2gxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i4 2083 d2gxdt(2,2,ilmn,ia,ispinor) = scale*buffer_r4 2084 d2gxdt(1,3,ilmn,ia,ispinor) =-scale*buffer_i5 2085 d2gxdt(2,3,ilmn,ia,ispinor) = scale*buffer_r5 2086 end if 2087 !$OMP END SINGLE 2088 else 2089 if (parity) then 2090 !$OMP SINGLE 2091 buffer_r1 = zero 2092 buffer_r2 = zero 2093 buffer_r3 = zero 2094 buffer_r4 = zero 2095 buffer_r5 = zero 2096 !$OMP END SINGLE 2097 !$OMP DO & 2098 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3,buffer_r4,buffer_r5), & 2099 !$OMP PRIVATE(aux_r) 2100 do ipw=1,npw 2101 aux_r = scalr(ipw)*ffnl(ipw,1+igamma,ilmn) 2102 buffer_r1 = buffer_r1 + aux_r 2103 buffer_r2 = buffer_r2 + aux_r*kpg(ipw,ibeta) 2104 aux_r = scalr(ipw)*ffnl(ipw,1+idelta,ilmn) 2105 buffer_r3 = buffer_r3 + aux_r 2106 buffer_r4 = buffer_r4 + aux_r*kpg(ipw,ibeta) 2107 aux_r = scalr(ipw)*ffnl(ipw,4+idelgam,ilmn) 2108 buffer_r5 = buffer_r5 + aux_r*kpg(ipw,ibeta) 2109 end do 2110 !$OMP SINGLE 2111 dgxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r1 2112 dgxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r3 2113 d2gxdt(1,1,ilmn,ia,ispinor) = scale*buffer_r2 2114 d2gxdt(1,2,ilmn,ia,ispinor) = scale*buffer_r4 2115 d2gxdt(1,3,ilmn,ia,ispinor) = scale*buffer_r5 2116 !$OMP END SINGLE 2117 else 2118 !$OMP SINGLE 2119 buffer_i1 = zero 2120 buffer_i2 = zero 2121 buffer_i3 = zero 2122 buffer_i4 = zero 2123 buffer_i5 = zero 2124 !$OMP END SINGLE 2125 !$OMP DO & 2126 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3,buffer_i4,buffer_i5), & 2127 !$OMP PRIVATE(aux_i) 2128 do ipw=1,npw 2129 aux_i = scali(ipw)*ffnl(ipw,1+igamma,ilmn) 2130 buffer_i1 = buffer_i1 + aux_i 2131 buffer_i2 = buffer_i2 + aux_i*kpg(ipw,ibeta) 2132 aux_i = scali(ipw)*ffnl(ipw,1+idelta,ilmn) 2133 buffer_i3 = buffer_i3 + aux_i 2134 buffer_i4 = buffer_i4 + aux_i*kpg(ipw,ibeta) 2135 aux_i = scali(ipw)*ffnl(ipw,4+idelgam,ilmn) 2136 buffer_i5 = buffer_i5 + aux_i*kpg(ipw,ibeta) 2137 end do 2138 !$OMP END DO 2139 !$OMP SINGLE 2140 dgxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i1 2141 dgxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i3 2142 d2gxdt(1,1,ilmn,ia,ispinor) =-scale*buffer_i2 2143 d2gxdt(1,2,ilmn,ia,ispinor) =-scale*buffer_i4 2144 d2gxdt(1,3,ilmn,ia,ispinor) =-scale*buffer_i5 2145 !$OMP END SINGLE 2146 end if 2147 end if 2148 end do 2149 end if 2150 2151 ! -------------------------------------------------------------------- 2152 ! CHOICE= 4, 24 -- SIGNS= 1 2153 ! Accumulate d2Gxdt --- 2nd derivative wrt 2 atm pos. --- for all dirs 2154 ! -------------------------------------------------------------------- 2155 if ((signs==1).and.(choice_==4.or.choice_==24)) then 2156 mu0=1 2157 do ilmn=1,nlmn 2158 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2159 scale=wt;if (il>1) scale=-scale 2160 scale2=two_pi2*scale 2161 if (cplex==2) then 2162 !$OMP SINGLE 2163 buffer_ra = zero ; buffer_rb = zero 2164 buffer_rc = zero ; buffer_rd = zero 2165 buffer_re = zero ; buffer_rf = zero 2166 buffer_ia = zero ; buffer_ib = zero 2167 buffer_ic = zero ; buffer_id = zero 2168 buffer_ie = zero ; buffer_if = zero 2169 !$OMP END SINGLE 2170 !$OMP DO & 2171 !$OMP PRIVATE(aux_r,aux_i), & 2172 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc,buffer_rd,buffer_re,buffer_rf),& 2173 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic,buffer_id,buffer_ie,buffer_if) 2174 do ipw=1,npw 2175 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 2176 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 2177 buffer_ra = buffer_ra - aux_r*kpg(ipw,4) 2178 buffer_rb = buffer_rb - aux_r*kpg(ipw,5) 2179 buffer_rc = buffer_rc - aux_r*kpg(ipw,6) 2180 buffer_rd = buffer_rd - aux_r*kpg(ipw,7) 2181 buffer_re = buffer_re - aux_r*kpg(ipw,8) 2182 buffer_rf = buffer_rf - aux_r*kpg(ipw,9) 2183 buffer_ia = buffer_ia - aux_i*kpg(ipw,4) 2184 buffer_ib = buffer_ib - aux_i*kpg(ipw,5) 2185 buffer_ic = buffer_ic - aux_i*kpg(ipw,6) 2186 buffer_id = buffer_id - aux_i*kpg(ipw,7) 2187 buffer_ie = buffer_ie - aux_i*kpg(ipw,8) 2188 buffer_if = buffer_if - aux_i*kpg(ipw,9) 2189 end do 2190 !$OMP END DO 2191 !$OMP SINGLE 2192 if (parity) then 2193 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ra 2194 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ia 2195 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb 2196 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_ib 2197 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc 2198 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_ic 2199 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd 2200 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_id 2201 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re 2202 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_ie 2203 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf 2204 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_if 2205 else 2206 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_ia 2207 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ra 2208 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib 2209 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb 2210 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic 2211 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc 2212 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id 2213 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd 2214 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie 2215 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re 2216 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if 2217 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf 2218 end if 2219 !$OMP END SINGLE 2220 else 2221 if (parity) then 2222 !$OMP SINGLE 2223 buffer_ra = zero ; buffer_rb = zero 2224 buffer_rc = zero ; buffer_rd = zero 2225 buffer_re = zero ; buffer_rf = zero 2226 !$OMP END SINGLE 2227 !$OMP DO & 2228 !$OMP PRIVATE(aux_r,aux_i), & 2229 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc,buffer_rd,buffer_re,buffer_rf) 2230 do ipw=1,npw 2231 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 2232 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 2233 buffer_ra = buffer_ra - aux_r*kpg(ipw,4) 2234 buffer_rb = buffer_rb - aux_r*kpg(ipw,5) 2235 buffer_rc = buffer_rc - aux_r*kpg(ipw,6) 2236 buffer_rd = buffer_rd - aux_r*kpg(ipw,7) 2237 buffer_re = buffer_re - aux_r*kpg(ipw,8) 2238 buffer_rf = buffer_rf - aux_r*kpg(ipw,9) 2239 end do 2240 !$OMP END DO 2241 !$OMP SINGLE 2242 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_ra 2243 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_rb 2244 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_rc 2245 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_rd 2246 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_re 2247 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_rf 2248 !$OMP END SINGLE 2249 else 2250 !$OMP SINGLE 2251 buffer_ia = zero ; buffer_ib = zero 2252 buffer_ic = zero ; buffer_id = zero 2253 buffer_ie = zero ; buffer_if = zero 2254 !$OMP END SINGLE 2255 !$OMP DO & 2256 !$OMP PRIVATE(aux_i,aux_r), & 2257 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic,buffer_id,buffer_ie,buffer_if) 2258 do ipw=1,npw 2259 aux_r = scalr(ipw)*ffnl(ipw,1,ilmn) 2260 aux_i = scali(ipw)*ffnl(ipw,1,ilmn) 2261 buffer_ia = buffer_ia - aux_i*kpg(ipw,4) 2262 buffer_ib = buffer_ib - aux_i*kpg(ipw,5) 2263 buffer_ic = buffer_ic - aux_i*kpg(ipw,6) 2264 buffer_id = buffer_id - aux_i*kpg(ipw,7) 2265 buffer_ie = buffer_ie - aux_i*kpg(ipw,8) 2266 buffer_if = buffer_if - aux_i*kpg(ipw,9) 2267 end do 2268 !$OMP END DO 2269 !$OMP SINGLE 2270 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_ia 2271 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_ib 2272 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_ic 2273 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_id 2274 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_ie 2275 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_if 2276 !$OMP END SINGLE 2277 end if 2278 end if 2279 end do 2280 end if 2281 2282 ! -------------------------------------------------------------------- 2283 ! CHOICE= 5, 51, 52, 53, 54, 55, 8, 81 -- SIGNS= 1 2284 ! Accumulate dGxdt --- derivative wrt k --- for all directions 2285 ! -------------------------------------------------------------------- 2286 if ((signs==1).and.& 2287 & (choice_==5 .or.choice_==51.or.choice_==52.or.choice_==53.or.& 2288 & choice_==54.or.choice_==55.or.choice_==8 .or.choice_==81)) then 2289 mu0=1 2290 if (choice_==54) mu0=4 2291 if (choice_==55) mu0=7 2292 do ilmn=1,nlmn 2293 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2294 scale=wt;if (il>1) scale=-scale 2295 if (cplex==2) then 2296 !$OMP SINGLE 2297 buffer_r1 = zero ; buffer_r2 = zero 2298 buffer_r3 = zero 2299 buffer_i1 = zero ; buffer_i2 = zero 2300 buffer_i3 = zero 2301 !$OMP END SINGLE 2302 !$OMP DO & 2303 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3), & 2304 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3) 2305 do ipw=1,npw 2306 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn) 2307 buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn) 2308 buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn) 2309 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn) 2310 buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn) 2311 buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn) 2312 end do 2313 !$OMP END DO 2314 !$OMP SINGLE 2315 if (parity) then 2316 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1 2317 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 2318 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2 2319 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 2320 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3 2321 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 2322 else 2323 dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_i1 2324 dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1 2325 dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2 2326 dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2 2327 dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3 2328 dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3 2329 end if 2330 !$OMP END SINGLE 2331 else 2332 cplex_dgxdt(mu0:mu0+2) = 2 ! Warning dgxdt is here pure imaginary 2333 if (parity) then 2334 !$OMP SINGLE 2335 buffer_i1 = zero ; buffer_i2 = zero 2336 buffer_i3 = zero 2337 !$OMP END SINGLE 2338 !$OMP DO & 2339 !$OMP REDUCTION(+:buffer_i1,buffer_i2,buffer_i3) 2340 do ipw=1,npw 2341 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn) 2342 buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn) 2343 buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn) 2344 end do 2345 !$OMP END DO 2346 !$OMP SINGLE 2347 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1 2348 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2 2349 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3 2350 !$OMP END SINGLE 2351 else 2352 !$OMP SINGLE 2353 buffer_r1 = zero ; buffer_r2 = zero 2354 buffer_r3 = zero 2355 !$OMP END SINGLE 2356 !$OMP DO & 2357 !$OMP REDUCTION(+:buffer_r1,buffer_r2,buffer_r3) 2358 do ipw=1,npw 2359 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn) 2360 buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn) 2361 buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn) 2362 end do 2363 !$OMP END DO 2364 !$OMP SINGLE 2365 dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1 2366 dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2 2367 dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3 2368 !$OMP END SINGLE 2369 end if 2370 end if 2371 end do 2372 end if 2373 2374 ! -------------------------------------------------------------------- 2375 ! CHOICE= 5, 51, 52, 53, 54, 8, 81 -- SIGNS= 2 2376 ! Accumulate dGxdt --- derivative wrt k --- 2377 ! for direction IDIR (CHOICE=5,51,52) 2378 ! for direction IDIR2 (CHOICE=54,81) 2379 ! for directions IDIR-1 & IDIR+1 (CHOICE=53) 2380 ! for 2 directions (CHOICE=8) 2381 ! -------------------------------------------------------------------- 2382 if ((signs==2).and.& 2383 & (choice_==5.or.choice_==51.or.choice_==52.or.choice_==53.or.choice_==54.or. & 2384 & choice_==8.or.choice_==81)) then 2385 mu0=1 2386 if (choice_==5.or.choice_==51.or.choice_==52) then 2387 ! We compute the derivative in one direction 2388 ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir 2389 ndir=1 2390 end if 2391 if (choice_==53) then 2392 ! Case 53 though need multiple directions and here ffnl will contain the 2393 ! derivative information in locations 2,3, and 4 corresponding to idir = 1, 2, 2394 ! and 3. Moreover, choice 53 needs the derivatives in direction idir+1 and idir-1. 2395 ! The parameter vector ffnl_dir_dat contains the necessary translations in 2396 ! locations 1,2 for idir=1; 3,4 for idir=2; and 5,6 for idir=3. 2397 ffnl_dir(1)=ffnl_dir_dat(2*idir-1);ffnl_dir(2)=ffnl_dir_dat(2*idir) 2398 ndir=2 2399 end if 2400 if (choice_==54) then 2401 ! We compute the derivative in one direction (the second one) 2402 ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir2(idir) 2403 ndir=1 2404 end if 2405 if (choice_==8) then 2406 ! We compute the derivative in two directions 2407 ! depending on the input idir argument 2408 ffnl_dir(1)=idir1(idir);ffnl_dir(2)=idir2(idir) 2409 ndir=2 2410 end if 2411 if (choice_==81) then 2412 ! We compute the derivative in one direction (the second one) 2413 ffnl_dir(1)=idir2(idir) 2414 ndir=1 2415 end if 2416 do mua=1,ndir 2417 mu=mu0+mua-1 2418 do ilmn=1,nlmn 2419 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2420 scale=wt;if (il>1) scale=-scale 2421 if (cplex==2) then 2422 !$OMP SINGLE 2423 buffer_r1 = zero ; buffer_i1 = zero 2424 !$OMP END SINGLE 2425 !$OMP DO & 2426 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 2427 do ipw=1,npw 2428 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 2429 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 2430 end do 2431 !$OMP END DO 2432 !$OMP SINGLE 2433 if (parity) then 2434 dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1 2435 dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_i1 2436 else 2437 dgxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1 2438 dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_r1 2439 end if 2440 !$OMP END SINGLE 2441 else 2442 !$OMP SINGLE 2443 cplex_dgxdt(mu) = 2 ! Warning dgxdt is here pure imaginary 2444 !$OMP END SINGLE 2445 if (parity) then 2446 !$OMP SINGLE 2447 buffer_i1 = zero 2448 !$OMP END SINGLE 2449 !$OMP DO & 2450 !$OMP REDUCTION(+:buffer_i1) 2451 do ipw=1,npw 2452 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 2453 end do 2454 !$OMP END DO 2455 !$OMP SINGLE 2456 dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1 2457 !$OMP END SINGLE 2458 else 2459 !$OMP SINGLE 2460 buffer_r1 = zero 2461 !$OMP END SINGLE 2462 !$OMP DO & 2463 !$OMP REDUCTION(+:buffer_r1) 2464 do ipw=1,npw 2465 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn) 2466 end do 2467 !$OMP SINGLE 2468 dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1 2469 !$OMP END SINGLE 2470 end if 2471 end if 2472 end do 2473 end do 2474 end if 2475 2476 ! -------------------------------------------------------------------- 2477 ! CHOICE= 54 -- SIGNS= 1 2478 ! Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- in all dirs 2479 ! -------------------------------------------------------------------- 2480 if ((signs==1).and.(choice_==54)) then 2481 mu0=1 2482 do ilmn=1,nlmn 2483 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2484 scale=wt;if (il>1) scale=-scale 2485 scale2=two_pi*scale 2486 if (cplex==2) then 2487 do mua=1,3 ! atm. pos 2488 do mub=1,3 ! k 2489 mu=mu0+(mub-1)+3*(mua-1) 2490 !$OMP SINGLE 2491 buffer_r1 = zero ; buffer_i1 = zero 2492 !$OMP END SINGLE 2493 !$OMP DO & 2494 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 2495 do ipw=1,npw 2496 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2497 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2498 end do 2499 !$OMP END DO 2500 !$OMP SINGLE 2501 if (parity) then 2502 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 2503 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2504 else 2505 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_r1 2506 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 2507 end if 2508 !$OMP END SINGLE 2509 end do 2510 end do 2511 else 2512 !$OMP SINGLE 2513 cplex_d2gxdt(mu0:mu0+8) = 2 ! Warning dgxdt is here pure imaginary 2514 !$OMP END SINGLE 2515 if (parity) then 2516 do mua=1,3 ! atm. pos 2517 do mub=1,3 ! k 2518 mu=mu0+(mub-1)+3*(mua-1) 2519 !$OMP SINGLE 2520 buffer_r1 = zero 2521 !$OMP END SINGLE 2522 !$OMP DO & 2523 !$OMP REDUCTION(+:buffer_r1) 2524 do ipw=1,npw 2525 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2526 end do 2527 !$OMP END DO 2528 !$OMP SINGLE 2529 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2530 !$OMP END SINGLE 2531 end do 2532 end do 2533 else 2534 do mua=1,3 ! atm. pos 2535 do mub=1,3 ! k 2536 mu=mu0+(mub-1)+3*(mua-1) 2537 !$OMP SINGLE 2538 buffer_i1 = zero 2539 !$OMP END SINGLE 2540 !$OMP DO & 2541 !$OMP REDUCTION(+:buffer_i1) 2542 do ipw=1,npw 2543 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2544 end do 2545 !$OMP END DO 2546 !$OMP SINGLE 2547 d2gxdt(1,mu,ilmn,ia,ispinor) = -scale2*buffer_i1 2548 !$OMP END SINGLE 2549 end do 2550 end do 2551 end if 2552 end if 2553 end do 2554 end if 2555 2556 ! -------------------------------------------------------------------- 2557 ! CHOICE= 54 -- SIGNS= 2 2558 ! Accumulate d2Gxdt --- 2nd derivative wrt k and atm. pos --- for direction IDIR 2559 ! -------------------------------------------------------------------- 2560 if ((signs==2).and.(choice_==54)) then 2561 mu0=1 2562 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 2563 mua=idir1(idir) ! atm. pos 2564 mub=1;if(dimffnl>2) mub=idir2(idir) ! k 2565 do ilmn=1,nlmn 2566 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2567 scale=wt;if (il>1) scale=-scale 2568 scale2=two_pi*scale 2569 if (cplex==2) then 2570 !$OMP SINGLE 2571 buffer_r1 = zero ; buffer_i1 = zero 2572 !$OMP END SINGLE 2573 !$OMP DO & 2574 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 2575 do ipw=1,npw 2576 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2577 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2578 end do 2579 !$OMP END DO 2580 !$OMP SINGLE 2581 if (parity) then 2582 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 2583 d2gxdt(2,mu0,ilmn,ia,ispinor) = scale2*buffer_r1 2584 else 2585 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale2*buffer_r1 2586 d2gxdt(2,mu0,ilmn,ia,ispinor) =-scale2*buffer_i1 2587 end if 2588 !$OMP END SINGLE 2589 else 2590 cplex_d2gxdt(mu0) = 2 ! Warning d2gxdt is here pure imaginary 2591 if (parity) then 2592 !$OMP SINGLE 2593 buffer_r1 = zero 2594 !$OMP END SINGLE 2595 !$OMP DO & 2596 !$OMP REDUCTION(+:buffer_r1) 2597 do ipw=1,npw 2598 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2599 end do 2600 !$OMP END DO 2601 !$OMP SINGLE 2602 d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 2603 !$OMP END SINGLE 2604 else 2605 !$OMP SINGLE 2606 buffer_i1 = zero 2607 !$OMP END SINGLE 2608 !$OMP DO & 2609 !$OMP REDUCTION(+:buffer_i1) 2610 do ipw=1,npw 2611 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+mub,ilmn)*kpg(ipw,mua) 2612 end do 2613 !$OMP SINGLE 2614 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 2615 !$OMP END SINGLE 2616 end if 2617 end if 2618 end do 2619 end if 2620 2621 ! -------------------------------------------------------------------- 2622 ! CHOICE= 55 -- SIGNS= 1 2623 ! Accumulate d2Gxdt --- 2nd derivative wrt k and strains --- for all dirs 2624 ! -------------------------------------------------------------------- 2625 if ((signs==1).and.(choice_==55)) then 2626 mu0=1 2627 do ilmn=1,nlmn 2628 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2629 scale=wt;if (il>1) scale=-scale 2630 if (cplex==2) then 2631 do mua=1,6 ! strain 2632 do mub=1,3 ! k 2633 mu=mu0+(mub-1)+3*(mua-1) 2634 !$OMP SINGLE 2635 buffer_r1 = zero ; buffer_i1 = zero 2636 !$OMP END SINGLE 2637 !$OMP DO & 2638 !$OMP REDUCTION(+:buffer_r1,buffer_i1),& 2639 !$OMP PRIVATE(aux_r) 2640 do ipw=1,npw 2641 aux_r=ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub) 2642 buffer_r1 = buffer_r1 + aux_r*scalr(ipw) 2643 buffer_i1 = buffer_i1 + aux_r*scali(ipw) 2644 end do 2645 !$OMP END DO 2646 !$OMP SINGLE 2647 if (parity) then 2648 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1 2649 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_i1 2650 else 2651 d2gxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_i1 2652 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale*buffer_r1 2653 end if 2654 !$OMP END SINGLE 2655 end do 2656 end do 2657 else 2658 !$OMP SINGLE 2659 cplex_d2gxdt(mu0:mu0+17) = 2 ! Warning d2gxdt is here pure imaginary 2660 !$OMP END SINGLE 2661 if (parity) then 2662 do mua=1,6 ! strain 2663 do mub=1,3 ! k 2664 mu=mu0+(mub-1)+3*(mua-1) 2665 !$OMP SINGLE 2666 buffer_i1 = zero 2667 !$OMP END SINGLE 2668 !$OMP DO & 2669 !$OMP REDUCTION(+:buffer_i1) 2670 do ipw=1,npw 2671 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub) 2672 end do 2673 !$OMP END DO 2674 !$OMP SINGLE 2675 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1 2676 !$OMP END SINGLE 2677 end do 2678 end do 2679 else 2680 do mua=1,6 ! strain 2681 do mub=1,3 ! k 2682 mu=mu0+(mub-1)+3*(mua-1) 2683 !$OMP SINGLE 2684 buffer_r1 = zero 2685 !$OMP END SINGLE 2686 !$OMP DO & 2687 !$OMP REDUCTION(+:buffer_r1) 2688 do ipw=1,npw 2689 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+mua ,ilmn)*kpg(ipw,mub) 2690 end do 2691 !$OMP END DO 2692 !$OMP SINGLE 2693 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_r1 2694 !$OMP END SINGLE 2695 end do 2696 end do 2697 end if 2698 end if 2699 end do 2700 end if 2701 2702 ! -------------------------------------------------------------------- 2703 ! CHOICE= 6 -- SIGNS= 1 2704 ! Accumulate d2Gxdt --- 2nd derivative wrt atm. pos and strain --- for all dirs 2705 ! -------------------------------------------------------------------- 2706 if ((signs==1).and.(choice_==6)) then 2707 mu0=1;if (choice_==6) mu0=37 2708 ABI_MALLOC(scalarr,(npw,4)) 2709 ABI_MALLOC(scalari,(npw,4)) 2710 do ilmn=1,nlmn 2711 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2712 scale=wt;if (il>1) scale=-scale 2713 scale2=pi*scale 2714 do mu=1,4 2715 !$OMP SINGLE 2716 do ipw=1,npw 2717 scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn) 2718 scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn) 2719 end do 2720 !$OMP END SINGLE 2721 end do 2722 if(cplex==2) then 2723 do mub=1,6 2724 nub1=alpha(mub);nub2=beta(mub) 2725 do mua=1,3 2726 mu=mu0+(mua-1)+3*(mub-1) 2727 buffer_r1 = zero ; buffer_i1 = zero 2728 !$OMP DO & 2729 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 2730 do ipw=1,npw 2731 buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + & 2732 & scalari(ipw,1+nub2)*kpg(ipw,nub1)) 2733 buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + & 2734 & scalarr(ipw,1+nub2)*kpg(ipw,nub1)) 2735 end do 2736 !$OMP END DO 2737 if (parity) then 2738 !$OMP SINGLE 2739 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1 2740 d2gxdt(2,mu,ilmn,ia,ispinor) =-scale2*buffer_r1 2741 !$OMP END SINGLE 2742 else 2743 !$OMP SINGLE 2744 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2745 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1 2746 !$OMP END SINGLE 2747 end if 2748 end do 2749 end do 2750 else 2751 if (parity) then 2752 do mub=1,6 2753 nub1=alpha(mub);nub2=beta(mub) 2754 do mua=1,3 2755 mu=mu0+(mua-1)+3*(mub-1) 2756 buffer_i1 = zero 2757 !$OMP DO & 2758 !$OMP REDUCTION(+:buffer_r1) 2759 do ipw=1,npw 2760 buffer_i1 = buffer_i1 + kpg(ipw,mua)*(scalari(ipw,1+nub1)*kpg(ipw,nub2) + & 2761 & scalari(ipw,1+nub2)*kpg(ipw,nub1)) 2762 end do 2763 !$OMP END DO 2764 !$OMP SINGLE 2765 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_i1 2766 !$OMP END SINGLE 2767 end do 2768 end do 2769 else 2770 do mub=1,6 2771 nub1=alpha(mub);nub2=beta(mub) 2772 do mua=1,3 2773 mu=mu0+(mua-1)+3*(mub-1) 2774 buffer_r1 = zero 2775 !$OMP DO & 2776 !$OMP REDUCTION(+:buffer_r1) 2777 do ipw=1,npw 2778 buffer_r1 = buffer_r1 + kpg(ipw,mua)*(scalarr(ipw,1+nub1)*kpg(ipw,nub2) + & 2779 & scalarr(ipw,1+nub2)*kpg(ipw,nub1)) 2780 end do 2781 !$OMP END DO 2782 !$OMP SINGLE 2783 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2784 !$OMP END SINGLE 2785 end do 2786 end do 2787 end if 2788 end if 2789 end do 2790 ABI_FREE(scalarr) 2791 ABI_FREE(scalari) 2792 end if 2793 2794 ! -------------------------------------------------------------------- 2795 ! CHOICE= 6 -- SIGNS= 1 2796 ! Accumulate d2Gxdt --- 2nd derivative wrt 2 strains --- for all dirs 2797 ! -------------------------------------------------------------------- 2798 if ((signs==1).and.(choice_==6)) then 2799 mu0=1 2800 ABI_MALLOC(scalarr,(npw,10)) 2801 ABI_MALLOC(scalari,(npw,10)) 2802 do ilmn=1,nlmn 2803 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2804 scale=wt;if (il>1) scale=-scale 2805 scale2=quarter*scale 2806 do mu=1,10 2807 !$OMP SINGLE 2808 do ipw=1,npw 2809 scalarr(ipw,mu)=scalr(ipw)*ffnl(ipw,mu,ilmn) 2810 scalari(ipw,mu)=scali(ipw)*ffnl(ipw,mu,ilmn) 2811 end do 2812 !$OMP END SINGLE 2813 end do 2814 if(cplex==2) then 2815 do mub=1,6 2816 nub1=alpha(mub);nub2=beta(mub) 2817 do mua=1,6 2818 mu=mu0+(mua-1)+6*(mub-1) 2819 nua1=alpha(mua);nua2=beta(mua) 2820 gama=gamma(nub1,nua1) 2821 gamb=gamma(nub2,nua1) 2822 gamc=gamma(nub1,nua2) 2823 gamd=gamma(nub2,nua2) 2824 buffer_r1 = zero ; buffer_i1 = zero 2825 !$OMP DO & 2826 !$OMP REDUCTION(+:buffer_r1,buffer_i1), & 2827 !$OMP PRIVATE(kpga,kpgb,kpgc,kpgd) 2828 do ipw=1,npw 2829 kpga=kpg(ipw,nub1) 2830 kpgb=kpg(ipw,nub2) 2831 kpgc=kpg(ipw,nua1) 2832 kpgd=kpg(ipw,nua2) 2833 buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd & 2834 & + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc 2835 buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd & 2836 & + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc 2837 end do 2838 !$OMP END DO 2839 if (parity) then 2840 !$OMP SINGLE 2841 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2842 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_i1 2843 !$OMP END SINGLE 2844 else 2845 !$OMP SINGLE 2846 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 2847 d2gxdt(2,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2848 !$OMP END SINGLE 2849 end if 2850 end do 2851 end do 2852 else 2853 if (parity) then 2854 do mub=1,6 2855 nub1=alpha(mub);nub2=beta(mub) 2856 do mua=1,6 2857 mu=mu0+(mua-1)+6*(mub-1) 2858 nua1=alpha(mua);nua2=beta(mua) 2859 gama=gamma(nub1,nua1) 2860 gamb=gamma(nub2,nua1) 2861 gamc=gamma(nub1,nua2) 2862 gamd=gamma(nub2,nua2) 2863 buffer_r1 = zero 2864 !$OMP DO & 2865 !$OMP REDUCTION(+:buffer_r1), & 2866 !$OMP PRIVATE(kpga,kpgb,kpgc,kpgd) 2867 do ipw=1,npw 2868 kpga=kpg(ipw,nub1) 2869 kpgb=kpg(ipw,nub2) 2870 kpgc=kpg(ipw,nua1) 2871 kpgd=kpg(ipw,nua2) 2872 buffer_r1 = buffer_r1 + scalarr(ipw,4+gama)*kpgb*kpgd+scalarr(ipw,4+gamb)*kpga*kpgd & 2873 & + scalarr(ipw,4+gamc)*kpgb*kpgc+scalarr(ipw,4+gamd)*kpga*kpgc 2874 end do 2875 !$OMP END DO 2876 !$OMP SINGLE 2877 d2gxdt(1,mu,ilmn,ia,ispinor) = scale2*buffer_r1 2878 !$OMP END SINGLE 2879 end do 2880 end do 2881 else 2882 do mub=1,6 2883 nub1=alpha(mub);nub2=beta(mub) 2884 do mua=1,6 2885 mu=mu0+(mua-1)+6*(mub-1) 2886 nua1=alpha(mua);nua2=beta(mua) 2887 gama=gamma(nub1,nua1) 2888 gamb=gamma(nub2,nua1) 2889 gamc=gamma(nub1,nua2) 2890 gamd=gamma(nub2,nua2) 2891 buffer_i1 = zero 2892 !$OMP DO & 2893 !$OMP REDUCTION(+:buffer_r1) 2894 do ipw=1,npw 2895 kpga=kpg(ipw,nub1) 2896 kpgb=kpg(ipw,nub2) 2897 kpgc=kpg(ipw,nua1) 2898 kpgd=kpg(ipw,nua2) 2899 buffer_i1 = buffer_i1 + scalari(ipw,4+gama)*kpgb*kpgd+scalari(ipw,4+gamb)*kpga*kpgd & 2900 & + scalari(ipw,4+gamc)*kpgb*kpgc+scalari(ipw,4+gamd)*kpga*kpgc 2901 end do 2902 !$OMP END DO 2903 !$OMP SINGLE 2904 d2gxdt(1,mu,ilmn,ia,ispinor) =-scale2*buffer_i1 2905 !$OMP END SINGLE 2906 end do 2907 end do 2908 end if 2909 end if 2910 end do 2911 ABI_FREE(scalarr) 2912 ABI_FREE(scalari) 2913 end if 2914 2915 ! -------------------------------------------------------------------- 2916 ! CHOICE= 8, 81 -- SIGNS= 1 2917 ! Accumulate d2Gxdt --- 2nd derivative wrt 2 k wave vect. --- for all dirs 2918 ! -------------------------------------------------------------------- 2919 if ((signs==1).and.(choice_==8.or.choice_==81)) then 2920 mu0=1 2921 do ilmn=1,nlmn 2922 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 2923 scale=wt;if (il>1) scale=-scale 2924 if (cplex==2) then 2925 !$OMP SINGLE 2926 buffer_ra = zero ; buffer_rb = zero 2927 buffer_rc = zero ; buffer_rd = zero 2928 buffer_re = zero ; buffer_rf = zero 2929 buffer_ia = zero ; buffer_ib = zero 2930 buffer_ic = zero ; buffer_id = zero 2931 buffer_ie = zero ; buffer_if = zero 2932 !$OMP END SINGLE 2933 !$OMP DO & 2934 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc), & 2935 !$OMP REDUCTION(+:buffer_rd,buffer_re,buffer_rf), & 2936 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic), & 2937 !$OMP REDUCTION(+:buffer_id,buffer_ie,buffer_if) 2938 do ipw=1,npw 2939 buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn) 2940 buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn) 2941 buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn) 2942 buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn) 2943 buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn) 2944 buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn) 2945 buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn) 2946 buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn) 2947 buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn) 2948 buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn) 2949 buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn) 2950 buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn) 2951 end do 2952 !$OMP END DO 2953 !$OMP SINGLE 2954 if (parity) then 2955 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_ra 2956 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_ia 2957 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb 2958 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_ib 2959 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc 2960 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_ic 2961 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd 2962 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_id 2963 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re 2964 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_ie 2965 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf 2966 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_if 2967 else 2968 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_ia 2969 d2gxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_ra 2970 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib 2971 d2gxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb 2972 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic 2973 d2gxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc 2974 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id 2975 d2gxdt(2,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd 2976 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie 2977 d2gxdt(2,mu0+4,ilmn,ia,ispinor) = scale*buffer_re 2978 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if 2979 d2gxdt(2,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf 2980 end if 2981 !$OMP END SINGLE 2982 else 2983 if (parity) then 2984 !$OMP SINGLE 2985 buffer_ra = zero ; buffer_rb = zero 2986 buffer_rc = zero ; buffer_rd = zero 2987 buffer_re = zero ; buffer_rf = zero 2988 !$OMP END SINGLE 2989 !$OMP DO & 2990 !$OMP REDUCTION(+:buffer_ra,buffer_rb,buffer_rc), & 2991 !$OMP REDUCTION(+:buffer_rd,buffer_re,buffer_rf) 2992 do ipw=1,npw 2993 buffer_ra = buffer_ra + scalr(ipw)*ffnl(ipw,5,ilmn) 2994 buffer_rb = buffer_rb + scalr(ipw)*ffnl(ipw,6,ilmn) 2995 buffer_rc = buffer_rc + scalr(ipw)*ffnl(ipw,7,ilmn) 2996 buffer_rd = buffer_rd + scalr(ipw)*ffnl(ipw,8,ilmn) 2997 buffer_re = buffer_re + scalr(ipw)*ffnl(ipw,9,ilmn) 2998 buffer_rf = buffer_rf + scalr(ipw)*ffnl(ipw,10,ilmn) 2999 end do 3000 !$OMP END DO 3001 !$OMP SINGLE 3002 d2gxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_ra 3003 d2gxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_rb 3004 d2gxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_rc 3005 d2gxdt(1,mu0+3,ilmn,ia,ispinor) = scale*buffer_rd 3006 d2gxdt(1,mu0+4,ilmn,ia,ispinor) = scale*buffer_re 3007 d2gxdt(1,mu0+5,ilmn,ia,ispinor) = scale*buffer_rf 3008 !$OMP END SINGLE 3009 else 3010 !$OMP SINGLE 3011 buffer_ia = zero ; buffer_ib = zero 3012 buffer_ic = zero ; buffer_id = zero 3013 buffer_ie = zero ; buffer_if = zero 3014 !$OMP END SINGLE 3015 !$OMP DO & 3016 !$OMP REDUCTION(+:buffer_ia,buffer_ib,buffer_ic), & 3017 !$OMP REDUCTION(+:buffer_id,buffer_ie,buffer_if) 3018 do ipw=1,npw 3019 buffer_ia = buffer_ia + scali(ipw)*ffnl(ipw,5,ilmn) 3020 buffer_ib = buffer_ib + scali(ipw)*ffnl(ipw,6,ilmn) 3021 buffer_ic = buffer_ic + scali(ipw)*ffnl(ipw,7,ilmn) 3022 buffer_id = buffer_id + scali(ipw)*ffnl(ipw,8,ilmn) 3023 buffer_ie = buffer_ie + scali(ipw)*ffnl(ipw,9,ilmn) 3024 buffer_if = buffer_if + scali(ipw)*ffnl(ipw,10,ilmn) 3025 end do 3026 !$OMP END DO 3027 !$OMP SINGLE 3028 d2gxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_ia 3029 d2gxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_ib 3030 d2gxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_ic 3031 d2gxdt(1,mu0+3,ilmn,ia,ispinor) =-scale*buffer_id 3032 d2gxdt(1,mu0+4,ilmn,ia,ispinor) =-scale*buffer_ie 3033 d2gxdt(1,mu0+5,ilmn,ia,ispinor) =-scale*buffer_if 3034 !$OMP END SINGLE 3035 end if 3036 end if 3037 end do 3038 end if 3039 3040 ! -------------------------------------------------------------------- 3041 ! CHOICE= 8, 81 -- SIGNS= 2 3042 ! Accumulate d2Gxdt --- 2nd derivative wrt k --- for direction IDIR (Voigt not.) 3043 ! -------------------------------------------------------------------- 3044 if ((signs==2).and.(choice_==8.or.choice_==81)) then 3045 mu0=1 3046 !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9) 3047 mua=idir1(idir);mub=idir2(idir) 3048 !idir->(Voigt) (11->1,22->2,33->3,32->4,31->5,21->6) 3049 ffnl_dir(1)=gamma(mua,mub) 3050 do ilmn=1,nlmn 3051 il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0) 3052 scale=wt;if (il>1) scale=-scale 3053 if (cplex==2) then 3054 !$OMP SINGLE 3055 buffer_r1 = zero ; buffer_i1 = zero 3056 !$OMP END SINGLE 3057 !$OMP DO & 3058 !$OMP REDUCTION(+:buffer_r1,buffer_i1) 3059 do ipw=1,npw 3060 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 3061 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 3062 end do 3063 !$OMP END DO 3064 !$OMP SINGLE 3065 if (parity) then 3066 d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 3067 d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_i1 3068 else 3069 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 3070 d2gxdt(2,mu0,ilmn,ia,ispinor) = scale*buffer_r1 3071 end if 3072 !$OMP END SINGLE 3073 else 3074 if (parity) then 3075 !$OMP SINGLE 3076 buffer_r1 = zero 3077 !$OMP END SINGLE 3078 !$OMP DO & 3079 !$OMP REDUCTION(+:buffer_r1) 3080 do ipw=1,npw 3081 buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 3082 end do 3083 !$OMP END DO 3084 !$OMP SINGLE 3085 d2gxdt(1,mu0,ilmn,ia,ispinor) = scale*buffer_r1 3086 !$OMP END SINGLE 3087 else 3088 !$OMP SINGLE 3089 buffer_i1 = zero 3090 !$OMP END SINGLE 3091 !$OMP DO & 3092 !$OMP REDUCTION(+:buffer_i1) 3093 do ipw=1,npw 3094 buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,4+ffnl_dir(1),ilmn) 3095 end do 3096 !$OMP SINGLE 3097 d2gxdt(1,mu0,ilmn,ia,ispinor) =-scale*buffer_i1 3098 !$OMP END SINGLE 3099 end if 3100 end if 3101 end do 3102 end if 3103 3104 ! -------------------------------------------------------------------- 3105 ! END CHOICES 3106 ! -------------------------------------------------------------------- 3107 3108 end do ! End loop on atoms 3109 3110 ! Deallocate temporary space 3111 !$OMP SECTIONS 3112 !$OMP SECTION 3113 ABI_FREE(scali) 3114 !$OMP SECTION 3115 ABI_FREE(scalr) 3116 !$OMP END SECTIONS 3117 3118 end do ! End loop on spinorial components 3119 !$OMP END PARALLEL 3120 3121 ! ========================================================================== 3122 end if 3123 3124 !Has to reduce arrays in case of FFT parallelization 3125 if (mpi_enreg%nproc_fft>1) then 3126 call timab(48,1,tsec) 3127 if (choice>=0) then 3128 call xmpi_sum(gx,mpi_enreg%comm_fft,ierr) 3129 end if 3130 if (choice_>1) then 3131 if (ndgxdt>0) then 3132 call xmpi_sum(dgxdt,mpi_enreg%comm_fft,ierr) 3133 end if 3134 if (nd2gxdt>0) then 3135 call xmpi_sum(d2gxdt,mpi_enreg%comm_fft,ierr) 3136 end if 3137 end if 3138 call timab(48,2,tsec) 3139 end if 3140 3141 end subroutine opernla_ylm