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