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