Branch data 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 : 61711416 : subroutine opernla_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnl,d2gxdt,dgxdt,ffnl,gx,&
148 : 30855708 : & ia3,idir,indlmn,istwf_k,kpg,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpg,nlmn,&
149 [ - + ]: 30855708 : & 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 : 30855708 : real(dp),allocatable :: scali(:),scalr(:),scalari(:,:),scalarr(:,:)
190 : : ! *************************************************************************
191 : :
192 [ - + ]: 30855708 : if (choice==-1) return
193 : :
194 : : !Useful variables
195 : 30855708 : choice_=abs(choice)
196 [ + + ]: 30855708 : wt=four_pi/sqrt(ucvol);if (cplex==1) wt=2.d0*wt
197 [ + + + + ]: 61244073 : ipw0=1;if (istwf_k==2.and.mpi_enreg%me_g0_fft==1) ipw0=2
198 [ + + + + : 38188329 : cplex_dgxdt(:) = 0 ; if (cplex == 1) cplex_dgxdt(:) = 1
+ + ]
199 [ + + + + : 32803089 : cplex_d2gxdt(:) = 0 ; if (cplex == 1) cplex_d2gxdt(:) = 1
- + ]
200 : 30855708 : nthreads=1
201 : :
202 : : !Check compatibility of options
203 [ + + ]: 30855708 : 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 [ + + - + ]: 3420903 : & choice_==6 .or.choice_==8 .or.choice_==81)
209 : 0 : ABI_CHECK(check,'BUG: choice not compatible (for signs=1)')
210 : : end if
211 [ + + ]: 30855708 : 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 [ + + + + : 27434805 : & choice_== 8.or.choice_==81)
+ + - + ]
217 : 0 : ABI_CHECK(check,'BUG: signs=2 not compatible with this choice')
218 : : end if
219 [ + + + + ]: 30855708 : 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 [ + + + - ]: 30579756 : 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 [ + + + - ]: 30330924 : 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 [ + + + + : 29957676 : 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 [ + + + + : 29848935 : 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 [ - + ]: 30855708 : 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 [ + + ]: 63913278 : do ispinor =1,nspinor
263 : 33057570 : ipwshft=(ispinor-1)*npw
264 : :
265 : : ! Allocate work space
266 [ - + + - : 99172710 : ABI_MALLOC(scali,(npw))
+ - + - -
+ - + ]
267 [ - + + - : 66115140 : ABI_MALLOC(scalr,(npw))
- + - + ]
268 : :
269 : : ! Loop on atoms (blocking)
270 [ + + ]: 79750756 : do ia=1,nincat
271 [ + + ]: 46693186 : iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
272 : : ! Compute c(g).exp(2pi.i.g.R)
273 [ + + ]: 16028041665 : do ipw=ipw0,npw
274 : 15981348479 : jpw=ipw+ipwshft
275 : 15981348479 : scalr(ipw)=(vect(1,jpw)*ph3d(1,ipw,iaph3d)-vect(2,jpw)*ph3d(2,ipw,iaph3d))
276 : 16028041665 : scali(ipw)=(vect(2,jpw)*ph3d(1,ipw,iaph3d)+vect(1,jpw)*ph3d(2,ipw,iaph3d))
277 : : end do
278 : :
279 [ + + ]: 46693186 : if (ipw0==2) then
280 : 1850060 : scalr(1)=half*vect(1,1+ipwshft)*ph3d(1,1,iaph3d)
281 : 1850060 : 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 [ + + ]: 46693186 : if (choice>=0) then
290 : :
291 [ + + ]: 548327010 : do ilmn=1,nlmn
292 : 501978231 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
293 [ + + ]: 501978231 : scale=wt;if (il>1) scale=-scale
294 [ + + ]: 548327010 : if (cplex==2) then
295 : : buffer_r1 = zero ; buffer_i1 = zero
296 [ + + ]: >12746*10^7 : do ipw=1,npw
297 : >12699*10^7 : buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)
298 : >12746*10^7 : buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)
299 : : end do
300 [ + + ]: 463072925 : if (parity) then
301 : 210398608 : gx(1,ilmn,ia,ispinor) = scale*buffer_r1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_i1
302 : : else
303 : 252674317 : gx(1,ilmn,ia,ispinor) =-scale*buffer_i1 ; gx(2,ilmn,ia,ispinor) = scale*buffer_r1
304 : : end if
305 : : else
306 [ + + ]: 38905306 : if (parity) then
307 : : buffer_r1 = zero
308 [ + + ]: 11186780386 : do ipw=1,npw
309 : 11186780386 : buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1,ilmn)
310 : : end do
311 : 11102704 : gx(1,ilmn,ia,ispinor) = scale*buffer_r1
312 : : else
313 : : buffer_i1 = zero
314 [ + + ]: 28712182242 : do ipw=1,npw
315 : 28712182242 : buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1,ilmn)
316 : : end do
317 : 27802602 : 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 [ + + + + ]: 46693186 : if ((signs==1).and. &
328 : : & (choice_==2.or.choice_==4.or.choice_==6.or.choice_==23.or.choice_==24.or.choice_==54)) then
329 [ + + ]: 1259538 : mu0=1;if (choice_==23.or.choice_==6) mu0=7
330 : :
331 [ + + ]: 13411421 : do ilmn=1,nlmn
332 : 12151883 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
333 [ + + ]: 12151883 : scale=wt;if (il>1) scale=-scale
334 : 12151883 : scale2=two_pi*scale
335 [ + + ]: 54761042 : 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 [ + + ]: 1964883120 : do ipw=1,npw
340 : 1956453373 : aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
341 : 1956453373 : aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
342 : 1956453373 : buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1)
343 : 1956453373 : buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2)
344 : 1956453373 : buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3)
345 : 1956453373 : buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1)
346 : 1956453373 : buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2)
347 : 1964883120 : buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3)
348 : : end do
349 [ + + ]: 8429747 : if (parity) then
350 : 3807271 : dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1
351 : 3807271 : dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_i1
352 : 3807271 : dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
353 : 3807271 : dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_i2
354 : 3807271 : dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
355 : 3807271 : dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_i3
356 : : else
357 : 4622476 : dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_i1
358 : 4622476 : dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1
359 : 4622476 : dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2
360 : 4622476 : dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
361 : 4622476 : dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale2*buffer_i3
362 : 4622476 : dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
363 : : end if
364 : : else
365 [ + + ]: 3722136 : if (parity) then
366 : : buffer_r1 = zero ; buffer_r2 = zero ; buffer_r3 = zero
367 [ + + ]: 1176050806 : do ipw=1,npw
368 : 1175077572 : aux_i = scali(ipw)*ffnl(ipw,1,ilmn)
369 : 1175077572 : buffer_r1 = buffer_r1 - aux_i*kpg(ipw,1)
370 : 1175077572 : buffer_r2 = buffer_r2 - aux_i*kpg(ipw,2)
371 : 1176050806 : buffer_r3 = buffer_r3 - aux_i*kpg(ipw,3)
372 : : end do
373 : 973234 : dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale2*buffer_r1
374 : 973234 : dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale2*buffer_r2
375 : 973234 : dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale2*buffer_r3
376 : : else
377 : : buffer_i1 = zero ; buffer_i2 = zero ; buffer_i3 = zero
378 [ + + ]: 3151478080 : do ipw=1,npw
379 : 3148729178 : aux_r = scalr(ipw)*ffnl(ipw,1,ilmn)
380 : 3148729178 : buffer_i1 = buffer_i1 + aux_r*kpg(ipw,1)
381 : 3148729178 : buffer_i2 = buffer_i2 + aux_r*kpg(ipw,2)
382 : 3151478080 : buffer_i3 = buffer_i3 + aux_r*kpg(ipw,3)
383 : : end do
384 : 2748902 : dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale2*buffer_i1
385 : 2748902 : dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale2*buffer_i2
386 : 2748902 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : if ((signs==1).and.(choice_==3.or.choice_==6.or.choice_==23.or.choice_==55)) then
488 : :
489 : 3076394 : mu0=1
490 [ + + ]: 3076394 : do ilmn=1,nlmn
491 : 2838077 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
492 [ + + ]: 2838077 : scale=wt;if (il>1) scale=-scale
493 : 2838077 : scale2=half*scale
494 [ + + ]: 44426015 : 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 [ + + ]: 1051849278 : do ipw=1,npw
502 : 1049076053 : aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn)
503 : 1049076053 : aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn)
504 : 1049076053 : aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn)
505 : 1049076053 : aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn)
506 : 1049076053 : aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn)
507 : 1049076053 : aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn)
508 : 1049076053 : buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1)
509 : 1049076053 : buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2)
510 : 1049076053 : buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3)
511 : 1049076053 : buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3)
512 : 1049076053 : buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3)
513 : 1049076053 : buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2)
514 : 1049076053 : buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1)
515 : 1049076053 : buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2)
516 : 1049076053 : buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3)
517 : 1049076053 : buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3)
518 : 1049076053 : buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3)
519 : 1051849278 : buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2)
520 : : end do
521 [ + + ]: 2773225 : if (parity) then
522 : 1359405 : dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1
523 : 1359405 : dgxdt(2,mu0 ,ilmn,ia,ispinor) =-scale*buffer_i1
524 : 1359405 : dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
525 : 1359405 : dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2
526 : 1359405 : dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
527 : 1359405 : dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3
528 : 1359405 : dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
529 : 1359405 : dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_i4
530 : 1359405 : dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
531 : 1359405 : dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_i5
532 : 1359405 : dgxdt(1,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
533 : 1359405 : dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_i6
534 : : else
535 : 1413820 : dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1
536 : 1413820 : dgxdt(2,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1
537 : 1413820 : dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
538 : 1413820 : dgxdt(2,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
539 : 1413820 : dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
540 : 1413820 : dgxdt(2,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
541 : 1413820 : dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4
542 : 1413820 : dgxdt(2,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
543 : 1413820 : dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5
544 : 1413820 : dgxdt(2,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
545 : 1413820 : dgxdt(1,mu0+5,ilmn,ia,ispinor) = scale2*buffer_i6
546 : 1413820 : dgxdt(2,mu0+5,ilmn,ia,ispinor) =-scale2*buffer_r6
547 : : end if
548 : : else
549 [ + + ]: 64852 : 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 [ + + ]: 37038778 : do ipw=1,npw
554 : 37007072 : aux_r2 = scalr(ipw)*ffnl(ipw,2,ilmn)
555 : 37007072 : aux_r3 = scalr(ipw)*ffnl(ipw,3,ilmn)
556 : 37007072 : aux_r4 = scalr(ipw)*ffnl(ipw,4,ilmn)
557 : 37007072 : buffer_r1 = buffer_r1 + aux_r2*kpg(ipw,1)
558 : 37007072 : buffer_r2 = buffer_r2 + aux_r3*kpg(ipw,2)
559 : 37007072 : buffer_r3 = buffer_r3 + aux_r4*kpg(ipw,3)
560 : 37007072 : buffer_r4 = buffer_r4 + aux_r4*kpg(ipw,2) + aux_r3*kpg(ipw,3)
561 : 37007072 : buffer_r5 = buffer_r5 + aux_r4*kpg(ipw,1) + aux_r2*kpg(ipw,3)
562 : 37038778 : buffer_r6 = buffer_r6 + aux_r3*kpg(ipw,1) + aux_r2*kpg(ipw,2)
563 : : end do
564 : 31706 : dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_r1
565 : 31706 : dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_r2
566 : 31706 : dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_r3
567 : 31706 : dgxdt(1,mu0+3,ilmn,ia,ispinor) =-scale2*buffer_r4
568 : 31706 : dgxdt(1,mu0+4,ilmn,ia,ispinor) =-scale2*buffer_r5
569 : 31706 : 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 [ + + ]: 34540264 : do ipw=1,npw
575 : 34507118 : aux_i2 = scali(ipw)*ffnl(ipw,2,ilmn)
576 : 34507118 : aux_i3 = scali(ipw)*ffnl(ipw,3,ilmn)
577 : 34507118 : aux_i4 = scali(ipw)*ffnl(ipw,4,ilmn)
578 : 34507118 : buffer_i1 = buffer_i1 + aux_i2*kpg(ipw,1)
579 : 34507118 : buffer_i2 = buffer_i2 + aux_i3*kpg(ipw,2)
580 : 34507118 : buffer_i3 = buffer_i3 + aux_i4*kpg(ipw,3)
581 : 34507118 : buffer_i4 = buffer_i4 + aux_i4*kpg(ipw,2) + aux_i3*kpg(ipw,3)
582 : 34507118 : buffer_i5 = buffer_i5 + aux_i4*kpg(ipw,1) + aux_i2*kpg(ipw,3)
583 : 34540264 : buffer_i6 = buffer_i6 + aux_i3*kpg(ipw,1) + aux_i2*kpg(ipw,2)
584 : : end do
585 : 33146 : dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1
586 : 33146 : dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
587 : 33146 : dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
588 : 33146 : dgxdt(1,mu0+3,ilmn,ia,ispinor) = scale2*buffer_i4
589 : 33146 : dgxdt(1,mu0+4,ilmn,ia,ispinor) = scale2*buffer_i5
590 : 33146 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + : 46693186 : 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 : 117616 : mu0=1
864 [ + + ]: 117616 : if (choice_==54) mu0=4
865 [ + + ]: 94192 : if (choice_==55) mu0=7
866 [ + + ]: 1076784 : do ilmn=1,nlmn
867 : 959168 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
868 [ + + ]: 959168 : scale=wt;if (il>1) scale=-scale
869 [ + - ]: 47652354 : 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 [ + + ]: 125910464 : do ipw=1,npw
875 : 124951296 : buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,2,ilmn)
876 : 124951296 : buffer_r2 = buffer_r2 + scalr(ipw)*ffnl(ipw,3,ilmn)
877 : 124951296 : buffer_r3 = buffer_r3 + scalr(ipw)*ffnl(ipw,4,ilmn)
878 : 124951296 : buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,2,ilmn)
879 : 124951296 : buffer_i2 = buffer_i2 + scali(ipw)*ffnl(ipw,3,ilmn)
880 : 125910464 : buffer_i3 = buffer_i3 + scali(ipw)*ffnl(ipw,4,ilmn)
881 : : end do
882 [ + + ]: 959168 : if (parity) then
883 : 250976 : dgxdt(1,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1
884 : 250976 : dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_i1
885 : 250976 : dgxdt(1,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
886 : 250976 : dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_i2
887 : 250976 : dgxdt(1,mu0+2,ilmn,ia,ispinor) = scale*buffer_r3
888 : 250976 : dgxdt(2,mu0+2,ilmn,ia,ispinor) = scale*buffer_i3
889 : : else
890 : 708192 : dgxdt(1,mu0 ,ilmn,ia,ispinor) =-scale*buffer_i1
891 : 708192 : dgxdt(2,mu0 ,ilmn,ia,ispinor) = scale*buffer_r1
892 : 708192 : dgxdt(1,mu0+1,ilmn,ia,ispinor) =-scale*buffer_i2
893 : 708192 : dgxdt(2,mu0+1,ilmn,ia,ispinor) = scale*buffer_r2
894 : 708192 : dgxdt(1,mu0+2,ilmn,ia,ispinor) =-scale*buffer_i3
895 : 708192 : 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 [ + + + + : 46693186 : 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 : 1703417 : mu0=1
938 [ + + - + ]: 1703417 : if (choice_==5.or.choice_==51.or.choice_==52) then
939 : : ! We compute the derivative in one direction
940 [ + + ]: 1573952 : ffnl_dir(1)=1;if(dimffnl>2) ffnl_dir(1)=idir
941 : : ndir=1
942 : : end if
943 [ - + ]: 1703417 : 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 [ + + ]: 1703417 : 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 [ + + ]: 1703417 : 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 [ + + ]: 1703417 : 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 [ + + ]: 3508651 : do mua=1,ndir
969 : 18272978 : mu=mu0+mua-1
970 [ + + ]: 64966164 : do ilmn=1,nlmn
971 : 16467744 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
972 [ + + ]: 16467744 : scale=wt;if (il>1) scale=-scale
973 [ + - ]: 18272978 : if (cplex==2) then
974 : : buffer_r1 = zero ; buffer_i1 = zero
975 [ + + ]: 2650987936 : do ipw=1,npw
976 : 2634520192 : buffer_r1 = buffer_r1 + scalr(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
977 : 2650987936 : buffer_i1 = buffer_i1 + scali(ipw)*ffnl(ipw,1+ffnl_dir(mua),ilmn)
978 : : end do
979 [ + + ]: 16467744 : if (parity) then
980 : 7658508 : dgxdt(1,mu,ilmn,ia,ispinor) = scale*buffer_r1
981 : 7658508 : dgxdt(2,mu,ilmn,ia,ispinor) = scale*buffer_i1
982 : : else
983 : 8809236 : dgxdt(1,mu,ilmn,ia,ispinor) =-scale*buffer_i1
984 : 8809236 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 46693186 : 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 [ + + + + ]: 79750756 : 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 [ - + ]: 33057570 : ABI_FREE(scali)
1489 [ - + ]: 63913278 : 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 [ + + ]: 30855708 : if (mpi_enreg%nproc_fft>1) then
3126 : 445964 : call timab(48,1,tsec)
3127 [ + + ]: 445964 : if (choice>=0) then
3128 : 443372 : call xmpi_sum(gx,mpi_enreg%comm_fft,ierr)
3129 : : end if
3130 [ + + ]: 445964 : 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 : 445964 : call timab(48,2,tsec)
3139 : : end if
3140 : :
3141 : : end subroutine opernla_ylm
3142 : : !!***
3143 : :
3144 : : end module m_opernla_ylm
3145 : : !!***
|
