Line data Source code
1 : !!****m* ABINIT/m_opernlb_ylm
2 : !! NAME
3 : !! m_opernlb_ylm
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-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_opernlb_ylm
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_errors
26 : #if defined HAVE_OPENMP
27 : use OMP_LIB
28 : #endif
29 :
30 : implicit none
31 :
32 : private
33 : !!***
34 :
35 : public :: opernlb_ylm
36 : integer,public,save :: opernlb_counter = -1
37 : !!***
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* ABINIT/opernlb_ylm
43 : !! NAME
44 : !! opernlb_ylm
45 : !!
46 : !! FUNCTION
47 : !! * Operate with the non-local part of the hamiltonian,
48 : !! from projected scalars to reciprocal space.
49 : !! * Operate with the non-local projectors and the overlap matrix,
50 : !! from projected scalars to reciprocal space.
51 : !!
52 : !! INPUTS
53 : !! choice=chooses possible output (see below)
54 : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
55 : !! 2 if <p_lmn|c> scalars are complex
56 : !! cplex_dgxdt(ndgxdt_fac) = used only when cplex = 1
57 : !! cplex_dgxdt(i)=1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
58 : !! cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
59 : !! dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Vnl (NL operator)
60 : !! dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfacrelated to Sij (overlap)
61 : !! dimffnl=second dimension of ffnl
62 : !! ffnl(npw,dimffnl,nlmn)= nonlocal quantities containing nonlocal form factors
63 : !! gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
64 : !! gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))= reduced projected scalars related to Sij (overlap)
65 : !! ia3=gives the number of the first atom in the subset presently treated
66 : !! idir=direction of the - atom to be moved in the case (choice=2,signs=2) or (choice=22,signs=2)
67 : !! - k point direction in the case (choice=5, 51, 52 and signs=2)
68 : !! - strain component (1:6) in the case (choice=2,signs=2) or (choice=6,signs=1)
69 : !! - strain component (1:9) in the case (choice=33,signs=2)
70 : !! - (1:9) components to specify the atom to be moved and the second q-gradient
71 : !! direction in the case (choice=25,signs=2)
72 : !! indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
73 : !! kpg(npw,nkpg)=(k+G) components (if nkpg=3).
74 : !! (k+G) Cartesian components for choice=33
75 : !! matblk=dimension of the array ph3d
76 : !! ndgxdtfac=second dimension of dgxdtfac
77 : !! nincat=number of atoms in the subset here treated
78 : !! nkpg=second dimension of array kpg (0 or 3)
79 : !! nlmn=number of (l,m,n) numbers for current type of atom
80 : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
81 : !! npw=number of plane waves in reciprocal space
82 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
83 : !! paw_opt= define the nonlocal operator concerned with:
84 : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
85 : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
86 : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
87 : !! paw_opt=3 : PAW overlap matrix (Sij)
88 : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
89 : !! ph3d(2,npw,matblk)=three-dimensional phase factors
90 : !! [qdir]= optional, direction of the q-gradient (only for choice=22, choice=25 and choice=33)
91 : !! ucvol=unit cell volume (bohr^3)
92 : !!
93 : !! OUTPUT
94 : !! (see side effects)
95 : !!
96 : !! SIDE EFFECTS
97 : !! --if (paw_opt=0)
98 : !! vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
99 : !! or one of its derivatives to the input vect.
100 : !! if (choice=22) <G|d2V_nonlocal/d(atm. pos)dq|vect_in> (at q=0)
101 : !! if (choice=25) <G|d3V_nonlocal/d(atm. pos)dqdq|vect_in> (at q=0)
102 : !! if (choice=33) <G|d2V_nonlocal/d(strain)dq|vect_in> (at q=0)
103 : !! --if (paw_opt=0, 1 or 4)
104 : !! vect(2,npwout*nspinor)=result of the aplication of the concerned operator
105 : !! or one of its derivatives to the input vect.:
106 : !! if (choice=1) <G|V_nonlocal|vect_in>
107 : !! if (choice=2) <G|dV_nonlocal/d(atm. pos)|vect_in>
108 : !! if (choice=3) <G|dV_nonlocal/d(strain)|vect_in>
109 : !! if (choice=5) <G|dV_nonlocal/d(k)|vect_in>
110 : !! if (choice=51) <G|d(right)V_nonlocal/d(k)|vect_in>
111 : !! if (choice=52) <G|d(left)V_nonlocal/d(k)|vect_in>
112 : !! if (choice=53) <G|d(twist)V_nonlocal/d(k)|vect_in>
113 : !! if (choice=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|vect_in>
114 : !! if (choice=8) <G|d2V_nonlocal/d(k)d(k)|vect_in>
115 : !! if (choice=81) <G|d[d(right)V_nonlocal/d(k)]/d(k)|vect_in>
116 : !! if (paw_opt=2)
117 : !! vect(2,npwout*nspinor)=final vector in reciprocal space:
118 : !! if (choice=1) <G|V_nonlocal-lamdba.(I+S)|vect_in> (note: not including <G|I|c>)
119 : !! if (choice=2) <G|d[V_nonlocal-lamdba.(I+S)]/d(atm. pos)|vect_in>
120 : !! if (choice=3) <G|d[V_nonlocal-lamdba.(I+S)]/d(strain)|vect_in>
121 : !! if (choice=5) <G|d[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
122 : !! if (choice=51) <G|d(right)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
123 : !! if (choice=52) <G|d(left)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
124 : !! if (choice=53) <G|d(twist)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
125 : !! if (choice=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|vect_in>
126 : !! if (choice=8) <G|d2[V_nonlocal-lamdba.(I+S)]/d(k)d(k)|vect_in>
127 : !! if (choice=81) <G|d[d(right[V_nonlocal-lamdba.(I+S)]/d(k)]/d(k)|vect_in>
128 : !! --if (paw_opt=3 or 4)
129 : !! svect(2,npwout*nspinor)=result of the aplication of Sij (overlap matrix)
130 : !! or one of its derivatives to the input vect.:
131 : !! if (choice=1) <G|I+S|vect_in> (note: not including <G|I|c>)
132 : !! if (choice=2) <G|dS/d(atm. pos)|vect_in>
133 : !! if (choice=3) <G|dS/d(strain)|vect_in>
134 : !! if (choice=5) <G|dS/d(k)|vect_in>
135 : !! if (choice=51) <G|d(right)S/d(k)|vect_in>
136 : !! if (choice=52) <G|d(left)S/d(k)|vect_in>
137 : !! if (choice=53) <G|d(twist)S/d(k)|vect_in>
138 : !! if (choice=54) <G|d[d(right)V_nonlocal/d(k)]/d(atm. pos)|vect_in>
139 : !! if (choice=7) <G|sum_i[p_i><p_i]|vect_in>
140 : !! if (choice=8) <G|d2S/d(k)d(k)|vect_in>
141 : !! if (choice=81) <G|d[d(right)S/d(k)]/d(k)|vect_in>
142 : !!
143 : !! NOTES
144 : !! 1-The openMP version is different from the standard version:
145 : !! the standard version is more effifient on one CPU core.
146 : !! 2-Operate for one type of atom, and within this given type of atom,
147 : !! for a subset of at most nincat atoms.
148 : !!
149 : !!
150 : !! SOURCE
151 :
152 35333413 : subroutine opernlb_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,&
153 35333413 : & d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnl,ffnl,gxfac,gxfac_sij,&
154 35333413 : & ia3,idir,indlmn,kpg,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpg,nlmn,nloalg,npw,&
155 35333413 : & nspinor,paw_opt,ph3d,svect,ucvol,vect,qdir)
156 :
157 : !Arguments ------------------------------------
158 : !scalars
159 : integer,intent(in) :: choice,cplex,cplex_fac,dimffnl,ia3,idir,matblk,ndgxdtfac,nd2gxdtfac,nincat
160 : integer,intent(in) :: nkpg,nlmn,npw,nspinor,paw_opt
161 : integer,intent(in),optional :: qdir
162 : real(dp),intent(in) :: ucvol
163 : !arrays
164 : integer,intent(in) :: cplex_dgxdt(ndgxdtfac),cplex_d2gxdt(nd2gxdtfac),indlmn(6,nlmn),nloalg(3)
165 : real(dp),intent(in) :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor)
166 : real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)
167 : real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat*(paw_opt/3),nspinor)
168 : real(dp),intent(in) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat*(paw_opt/3),nspinor)
169 : real(dp),intent(in) :: ffnl(npw,dimffnl,nlmn),gxfac(cplex_fac,nlmn,nincat,nspinor)
170 : real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
171 : real(dp),intent(in) :: kpg(npw,nkpg),ph3d(2,npw,matblk)
172 : real(dp),intent(inout) :: svect(:,:),vect(:,:)
173 : !Local variables-------------------------------
174 : !Arrays
175 : !scalars
176 : integer :: ia,ialpha,iaph3d,ibeta,ic,idelta,idelgam,igamma
177 : integer :: ii,il,ilmn,ipw,ipwshft,ispinor,jc,nthreads,ffnl_dir1,ffnl_dir(3)
178 : real(dp) :: scale,two_piinv,wt
179 : logical :: parity
180 : !arrays
181 : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
182 : integer,parameter :: idir1(9)=(/1,1,1,2,2,2,3,3,3/),idir2(9)=(/1,2,3,1,2,3,1,2,3/)
183 : integer,parameter :: nalpha(9)=(/1,2,3,3,3,2,2,1,1/),nbeta(9)=(/1,2,3,2,1,1,3,3,2/)
184 35333413 : real(dp),allocatable :: d2gxdtfac_(:,:,:),d2gxdtfacs_(:,:,:),dgxdtfac_(:,:,:),dgxdtfacs_(:,:,:),gxfac_(:,:),gxfacs_(:,:)
185 : ! real(dp),allocatable :: kpg(:,:)
186 35333413 : complex(dp),allocatable :: ztab(:)
187 : ! *************************************************************************
188 :
189 : DBG_ENTER("COLL")
190 :
191 : !Nothing to do when choice=4, 6 or 23
192 35333413 : if (choice==4.or.choice==6.or.choice==23) return
193 :
194 : !DDK not compatible with istwkf > 1
195 76368126 : if(cplex==1.and.(any(cplex_dgxdt(:)==2).or.any(cplex_d2gxdt(:)==2)))then
196 0 : ABI_BUG("opernlb_ylm+ddk not compatible with istwfk>1")
197 : end if
198 :
199 : !Inits
200 35333413 : wt=four_pi/sqrt(ucvol)
201 35333413 : nthreads=1
202 : #if defined HAVE_OPENMP
203 : nthreads=OMP_GET_NUM_THREADS()
204 : #endif
205 :
206 35333413 : if (paw_opt/=3) then
207 95979435 : ABI_MALLOC(gxfac_,(2,nlmn))
208 1040227496 : gxfac_(:,:)=zero
209 31993145 : if (choice>1) then
210 15773048 : ABI_MALLOC(dgxdtfac_,(2,ndgxdtfac,nlmn))
211 161042594 : if(ndgxdtfac>0) dgxdtfac_(:,:,:)=zero
212 : end if
213 31993145 : if (choice==54.or.choice==8.or.choice==81.or.choice==33) then
214 1812756 : ABI_MALLOC(d2gxdtfac_,(2,nd2gxdtfac,nlmn))
215 35711877 : if(nd2gxdtfac>0) d2gxdtfac_(:,:,:)=zero
216 : end if
217 : end if
218 35333413 : if (paw_opt>=3) then
219 54556215 : ABI_MALLOC(gxfacs_,(2,nlmn))
220 618520316 : gxfacs_(:,:)=zero
221 18185405 : if (choice>1) then
222 6800244 : ABI_MALLOC(dgxdtfacs_,(2,ndgxdtfac,nlmn))
223 61600860 : if (ndgxdtfac>0) dgxdtfacs_(:,:,:)=zero
224 : end if
225 18185405 : if (choice==54.or.choice==8.or.choice==81) then
226 131892 : ABI_MALLOC(d2gxdtfacs_,(2,nd2gxdtfac,nlmn))
227 1094769 : if (nd2gxdtfac>0) d2gxdtfacs_(:,:,:)=zero
228 : end if
229 : end if
230 :
231 : if (choice==33) two_piinv=1.0_dp/two_pi
232 :
233 35333413 : if (opernlb_counter>=0) then
234 0 : opernlb_counter = opernlb_counter + 1
235 0 : if (paw_opt==4) opernlb_counter = opernlb_counter + 1
236 : end if
237 :
238 106000239 : ABI_MALLOC(ztab,(npw))
239 :
240 : !==========================================================================
241 : !========== STANDARD VERSION ==============================================
242 : !==========================================================================
243 : if (nthreads==1) then
244 :
245 : ! Loop on spinorial components
246 72604850 : do ispinor=1,nspinor
247 37271437 : ipwshft=(ispinor-1)*npw
248 :
249 : ! Loop on atoms (blocking)
250 122456673 : do ia=1,nincat
251 49851823 : iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
252 : ! Scale gxfac with 4pi/sqr(omega).(-i)^l
253 49851823 : if (paw_opt/=3) then
254 521355787 : do ilmn=1,nlmn
255 476215948 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
256 476215948 : scale=wt;if (il>1) scale=-scale
257 521355787 : if (parity) then
258 601054463 : gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor)
259 203358273 : if (cplex_fac==1) gxfac_(2,ilmn)=zero
260 : else
261 272857675 : gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor)
262 272857675 : if (cplex_fac==2) then
263 250466579 : gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor)
264 : else
265 22391096 : gxfac_(1,ilmn)=zero
266 : end if
267 : end if! parity
268 : end do ! ilmn
269 45139839 : if (choice>1) then
270 51673360 : do ilmn=1,nlmn
271 46592643 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
272 46592643 : scale=wt;if (il>1) scale=-scale
273 51673360 : if (parity) then
274 20173557 : if(cplex_fac==2)then
275 92582031 : dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor)
276 : else
277 0 : do ii=1,ndgxdtfac
278 0 : ic = cplex_dgxdt(ii) ; jc = 3-ic
279 0 : dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
280 0 : dgxdtfac_(jc,ii,ilmn)=zero
281 : end do
282 : end if
283 : else
284 26419086 : if(cplex_fac==2)then
285 56311506 : do ii=1,ndgxdtfac
286 29892420 : dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor)
287 56311506 : dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
288 : end do
289 : else
290 0 : do ii=1,ndgxdtfac
291 0 : ic = cplex_dgxdt(ii) ; jc = 3-ic
292 0 : dgxdtfac_(ic,ii,ilmn)=zero
293 0 : if(ic==1)then
294 0 : dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
295 : else
296 0 : dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
297 : end if
298 : end do
299 : end if
300 : end if
301 : end do
302 : end if ! choice>1
303 45139839 : if (choice==54.or.choice==8.or.choice==81.or.choice==33) then
304 8285016 : do ilmn=1,nlmn
305 7435935 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
306 7435935 : scale=wt;if (il>1) scale=-scale
307 8285016 : if (parity) then
308 3962601 : if(cplex_fac==2)then
309 36752292 : d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor)
310 : else
311 0 : do ii=1,nd2gxdtfac
312 0 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
313 0 : d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
314 0 : d2gxdtfac_(jc,ii,ilmn)=zero
315 : end do
316 : end if
317 : else
318 3473334 : if(cplex_fac==2)then
319 12918636 : do ii=1,nd2gxdtfac
320 9445302 : d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor)
321 12918636 : d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
322 : end do
323 : else
324 0 : do ii=1,nd2gxdtfac
325 0 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
326 0 : d2gxdtfac_(ic,ii,ilmn)=zero
327 0 : if(ic==1)then
328 0 : d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
329 : else
330 0 : d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
331 : end if
332 : end do
333 : end if
334 : end if
335 : end do! ilmn
336 : end if ! choice 54 or 8 or 81 or 33
337 : end if ! paw_opt /= 3
338 :
339 : ! Scale gxfac_sij with 4pi/sqr(omega).(-i)^l
340 49851823 : if (paw_opt>=3) then
341 :
342 342576453 : do ilmn=1,nlmn
343 313789454 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
344 313789454 : scale=wt;if (il>1) scale=-scale
345 342576453 : if (parity) then
346 404164864 : gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor)
347 137698424 : if (cplex==1) gxfacs_(2,ilmn)=zero
348 : else
349 176091030 : gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor)
350 176091030 : if (cplex==2) then
351 153252314 : gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor)
352 : else
353 22838716 : gxfacs_(1,ilmn)=zero
354 : end if
355 : end if
356 : end do
357 28786999 : if (choice>1) then
358 18985844 : do ilmn=1,nlmn
359 17036738 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
360 17036738 : scale=wt;if (il>1) scale=-scale
361 18985844 : if (parity) then
362 5281730 : if(cplex==2)then
363 21175286 : dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn) = scale * dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor)
364 : else
365 0 : do ii=1,ndgxdtfac
366 0 : ic = cplex_dgxdt(ii) ; jc = 3-ic
367 0 : dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
368 0 : dgxdtfacs_(jc,ii,ilmn)=zero
369 : end do
370 : end if
371 : else
372 11755008 : if(cplex==2)then
373 23540292 : do ii=1,ndgxdtfac
374 11785284 : dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor)
375 23540292 : dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
376 : end do
377 : else
378 0 : do ii=1,ndgxdtfac
379 0 : ic = cplex_dgxdt(ii) ; jc = 3-ic
380 0 : dgxdtfacs_(ic,ii,ilmn)=zero
381 0 : if(ic==1)then
382 0 : dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
383 : else
384 0 : dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
385 : end if
386 : end do
387 : end if
388 : end if
389 : end do ! ilmn
390 : end if ! choice>1
391 28786999 : if (choice==54.or.choice==8.or.choice==81) then
392 310644 : do ilmn=1,nlmn
393 276798 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
394 276798 : scale=wt;if (il>1) scale=-scale
395 29063797 : if (parity) then
396 73722 : if(cplex==2)then
397 294888 : d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor)
398 : else
399 0 : do ii=1,nd2gxdtfac
400 0 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
401 0 : d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
402 0 : d2gxdtfacs_(jc,ii,ilmn)=zero
403 : end do
404 : end if
405 : else
406 203076 : if(cplex==2)then
407 406152 : do ii=1,nd2gxdtfac
408 203076 : d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor)
409 406152 : d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
410 : end do
411 : else
412 0 : do ii=1,nd2gxdtfac
413 0 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
414 0 : d2gxdtfacs_(ic,ii,ilmn)=zero
415 0 : if(ic==1)then
416 0 : d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
417 : else
418 0 : d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
419 : end if
420 : end do
421 : end if
422 : end if ! parity
423 : end do ! ilmn
424 : end if ! choice == 54 or 8 or 81
425 : end if ! paw_opt >= 3
426 :
427 : ! Compute <g|Vnl|c> (or derivatives) for each plane wave:
428 :
429 49851823 : if (paw_opt/=3) then
430 :
431 13989613660 : ztab(:)=czero
432 :
433 : ! ------
434 45139839 : if (choice==1) then ! <g|Vnl|c>
435 469682427 : do ilmn=1,nlmn
436 >13744*10^7 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
437 : end do
438 : end if
439 :
440 : ! ------
441 45139839 : if (choice==2) then ! derivative w.r.t. atm. pos
442 14863658 : do ilmn=1,nlmn
443 1959612387 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp)
444 : end do
445 211704505 : ztab(:)=two_pi*kpg(:,idir)*ztab(:)
446 14863658 : do ilmn=1,nlmn
447 1959612387 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
448 : end do
449 : end if
450 :
451 : ! ------
452 45139839 : if (choice==22) then ! mixed derivative w.r.t. atm. pos and q vector (at q=0)
453 380000 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+qdir
454 380000 : if (idir==qdir) then
455 1279328 : do ilmn=1,nlmn
456 118987768 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
457 : end do
458 : end if
459 3647456 : do ilmn=1,nlmn
460 334731708 : ztab(:)=ztab(:)+kpg(:,idir)*ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
461 335111708 : ztab(:)=ztab(:)-ffnl(:,ffnl_dir1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
462 : end do
463 39224044 : ztab(:)=ztab(:)*two_pi
464 : end if
465 :
466 : ! ------
467 45139839 : if (choice==25) then ! mixed derivative w.r.t. atm. pos and two q vectors (at q=0)
468 : !Use same notation as the notes for clarity
469 248832 : ialpha=nalpha(idir)
470 248832 : idelta=nbeta(idir)
471 248832 : igamma=qdir
472 248832 : idelgam=gamma(idelta,igamma)
473 248832 : if (ialpha==igamma) then
474 801792 : do ilmn=1,nlmn
475 71810496 : ztab(:)=ztab(:)+ffnl(:,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
476 : end do
477 : end if
478 248832 : if (ialpha==idelta) then
479 801792 : do ilmn=1,nlmn
480 71810496 : ztab(:)=ztab(:)+ffnl(:,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
481 : end do
482 : end if
483 2405376 : do ilmn=1,nlmn
484 215182656 : ztab(:)=ztab(:)+kpg(:,ialpha)*ffnl(:,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
485 215431488 : ztab(:)=ztab(:)-ffnl(:,4+idelgam,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
486 : end do
487 24828768 : ztab(:)=ztab(:)*two_pi
488 : end if
489 :
490 : ! ------
491 45139839 : if (choice==3) then ! derivative w.r.t. strain
492 505584 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
493 505584 : if (idir<=3) then
494 2401476 : do ilmn=1,nlmn
495 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)&
496 : & *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp)&
497 447159214 : & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
498 : end do
499 : else
500 2427576 : do ilmn=1,nlmn
501 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)&
502 456186658 : & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
503 : end do
504 : end if
505 : end if
506 :
507 : ! ------
508 45139839 : if (choice==33) then ! mixed derivative w.r.t. strain and q vector (at q=0)
509 : !Use same notation as the notes for clarity
510 746496 : ibeta=nalpha(idir)
511 746496 : idelta=nbeta(idir)
512 746496 : igamma=qdir
513 746496 : idelgam=gamma(idelta,igamma)
514 746496 : if (ibeta==igamma) then
515 2405376 : do ilmn=1,nlmn
516 215182656 : ztab(:)=ztab(:)+onehalf*ffnl(:,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
517 215431488 : ztab(:)=ztab(:)+half*ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp)
518 : end do
519 : end if
520 746496 : if (ibeta==idelta) then
521 2405376 : do ilmn=1,nlmn
522 215182656 : ztab(:)=ztab(:)+onehalf*ffnl(:,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
523 215431488 : ztab(:)=ztab(:)+half*ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
524 : end do
525 : end if
526 7216128 : do ilmn=1,nlmn
527 645547968 : ztab(:)=ztab(:)+kpg(:,ibeta)*ffnl(:,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
528 645547968 : ztab(:)=ztab(:)+ffnl(:,1+idelta,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)
529 645547968 : ztab(:)=ztab(:)+ffnl(:,1+igamma,ilmn)*cmplx(d2gxdtfac_(1,2,ilmn),d2gxdtfac_(2,2,ilmn),kind=dp)
530 646294464 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfac_(1,3,ilmn),d2gxdtfac_(2,3,ilmn),kind=dp)
531 : end do
532 74486304 : ztab(:)=ztab(:)*two_piinv
533 : end if
534 :
535 : ! ------
536 45139839 : if (choice==5) then ! full derivative w.r.t. k
537 1701424 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
538 17642802 : do ilmn=1,nlmn
539 : ztab(:)=ztab(:)+ffnl(:,1 ,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)&
540 2550004980 : & +ffnl(:,ffnl_dir1,ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2 ,ilmn),kind=dp)
541 : end do
542 : end if
543 :
544 : ! ------
545 45139839 : if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi>
546 0 : do ilmn=1,nlmn
547 0 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
548 : end do
549 : end if
550 :
551 : ! ------
552 45139839 : if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi>
553 0 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
554 0 : do ilmn=1,nlmn
555 0 : ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
556 : end do
557 : end if
558 :
559 : ! ------
560 45139839 : if (choice==53) then ! twist derivative: <G|dp_i/dk_(idir+1)>V_ij<dp_j/dk_(idir+2)|psi>
561 0 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=(mod(idir,3)+1)+1
562 0 : do ilmn=1,nlmn
563 0 : il=indlmn(1,ilmn)
564 : ztab(:)=ztab(:) + &
565 0 : & ffnl(:,ffnl_dir1,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp)*CONJG(j_dpc**il)*four_pi
566 : end do
567 : end if
568 :
569 : ! ------
570 45139839 : if (choice==54) then ! mixed derivative w.r.t. atm. pos and (right) k
571 0 : do ilmn=1,nlmn
572 0 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp)
573 : end do
574 0 : ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:)
575 0 : do ilmn=1,nlmn
576 0 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)
577 : end do
578 : end if
579 :
580 : ! ------
581 45139839 : if (choice==8) then ! full second order derivative w.r.t. k
582 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
583 102585 : ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
584 102585 : ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
585 1068888 : do ilmn=1,nlmn
586 : ztab(:)=ztab(:) &
587 : & +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2, ilmn),kind=dp)&
588 : & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
589 : & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)&
590 163680645 : & +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
591 : end do
592 : end if
593 :
594 : ! ------
595 45139839 : if (choice==81) then
596 : ! partial second order derivative w.r.t. k
597 : ! full derivative w.r.t. k1, right derivative w.r.t. k2
598 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
599 0 : ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
600 0 : do ilmn=1,nlmn
601 : ztab(:)=ztab(:) &
602 : & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
603 0 : & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
604 : end do
605 : end if
606 :
607 : ! ------
608 13989613660 : ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
609 :
610 13989613660 : vect(1,1+ipwshft:npw+ipwshft)=vect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
611 13989613660 : vect(2,1+ipwshft:npw+ipwshft)=vect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:))
612 :
613 : end if ! paw_opt /= 3
614 :
615 : ! Compute <g|S|c> (or derivatives) for each plane wave:
616 :
617 87123260 : if (paw_opt>=3) then
618 :
619 11519367342 : ztab(:)=czero
620 :
621 : ! ------
622 28786999 : if (choice==1) then ! <g|S|c>
623 323590609 : do ilmn=1,nlmn
624 >11744*10^7 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
625 : end do
626 : end if
627 :
628 : ! ------
629 28786999 : if (choice==2) then ! derivative w.r.t. atm. pos
630 7173664 : do ilmn=1,nlmn
631 827164420 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp)
632 : end do
633 95265468 : ztab(:)=two_pi*kpg(:,idir)*ztab(:)
634 7173664 : do ilmn=1,nlmn
635 827164420 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
636 : end do
637 : end if
638 :
639 : ! ------
640 28786999 : if (choice==3) then ! derivative w.r.t. strain
641 304944 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
642 304944 : if (idir<=3) then
643 1409172 : do ilmn=1,nlmn
644 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)&
645 : & *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)&
646 371926036 : & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
647 : end do
648 : else
649 1435164 : do ilmn=1,nlmn
650 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)&
651 380941468 : & -ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
652 : end do
653 : end if
654 : end if
655 :
656 : ! ------
657 28786999 : if (choice==5) then ! full derivative w.r.t. k
658 613832 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
659 6387208 : do ilmn=1,nlmn
660 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)&
661 952530360 : & +ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
662 : end do
663 : end if
664 :
665 : ! ------
666 28786999 : if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi>
667 2269992 : do ilmn=1,nlmn
668 379979152 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
669 : end do
670 : end if
671 :
672 : ! ------
673 28786999 : if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi>
674 0 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
675 0 : do ilmn=1,nlmn
676 0 : ztab(:)=ztab(:)+ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
677 : end do
678 : end if
679 :
680 : ! ------
681 28786999 : if (choice==53) then ! twist derivative: <G|dp_i/dk_(idir+1)>S_ij<dp_j/dk_(idir+2)|psi>
682 0 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=(mod(idir,3)+1)+1
683 0 : do ilmn=1,nlmn
684 0 : il=indlmn(1,ilmn)
685 : ztab(:)=ztab(:) + &
686 0 : & ffnl(:,ffnl_dir1,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp)*CONJG(j_dpc**il)*four_pi
687 : end do
688 : end if
689 :
690 : ! ------
691 28786999 : if (choice==54) then ! mixed derivative w.r.t. atm. pos and k
692 207360 : do ilmn=1,nlmn
693 17579520 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp)
694 : end do
695 2194560 : ztab(:)=two_pi*kpg(:,idir1(idir))*ztab(:)
696 207360 : do ilmn=1,nlmn
697 17579520 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)
698 : end do
699 : end if
700 :
701 : ! ------
702 28786999 : if (choice==8) then ! full second order derivative w.r.t. k
703 5046 : ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
704 5046 : ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
705 51444 : do ilmn=1,nlmn
706 : ztab(:)=ztab(:) &
707 : & +ffnl(:,4+ffnl_dir(3),ilmn)*cmplx( gxfacs_(1, ilmn), gxfacs_(2, ilmn),kind=dp)&
708 : & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
709 : & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)&
710 5125110 : & +ffnl(:,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
711 : end do
712 : end if
713 :
714 : ! ------
715 28786999 : if (choice==81) then
716 : ! partial second order derivative w.r.t. k
717 : ! full derivative w.r.t. k1, right derivative w.r.t. k2
718 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
719 5760 : ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
720 51840 : do ilmn=1,nlmn
721 : ztab(:)=ztab(:) &
722 : & +ffnl(:,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
723 4394880 : & +ffnl(:,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
724 : end do
725 : end if
726 :
727 :
728 : ! ------
729 11519367342 : ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
730 11519367342 : svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
731 11519367342 : svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:))
732 : end if ! paw_opt >= 3
733 :
734 : ! End loop on atoms
735 : end do
736 : end do ! End loop on spinors
737 :
738 :
739 : ! ==========================================================================
740 : ! ========== OPENMP VERSION ================================================
741 : ! ==========================================================================
742 : else
743 :
744 : ! Loop on spinorial components
745 : do ispinor=1,nspinor
746 : ipwshft=(ispinor-1)*npw
747 :
748 : ! Loop on atoms (blocking)
749 : do ia=1,nincat
750 : iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
751 :
752 : ! Scale gxfac with 4pi/sqr(omega).(-i)^l
753 : if (paw_opt/=3) then
754 : !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc)
755 : !$OMP DO
756 : do ilmn=1,nlmn
757 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
758 : scale=wt;if (il>1) scale=-scale
759 : if (parity) then
760 : gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor)
761 : if (cplex_fac==1) gxfac_(2,ilmn)=zero
762 : else
763 : gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor)
764 : if (cplex_fac==2) then
765 : gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor)
766 : else
767 : gxfac_(1,ilmn)=zero
768 : end if
769 : end if
770 : end do
771 : !$OMP END DO
772 : if (choice>1) then
773 : !$OMP DO
774 : do ilmn=1,nlmn
775 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
776 : scale=wt;if (il>1) scale=-scale
777 : if (parity) then
778 : if(cplex_fac==2)then
779 : dgxdtfac_(1:cplex_fac,1:ndgxdtfac,ilmn)=scale*dgxdtfac(1:cplex_fac,1:ndgxdtfac,ilmn,ia,ispinor)
780 : else
781 : do ii=1,ndgxdtfac
782 : ic = cplex_dgxdt(ii) ; jc = 3-ic
783 : dgxdtfac_(ic,ii,ilmn)=scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
784 : dgxdtfac_(jc,ii,ilmn)=zero
785 : end do
786 : end if
787 : else
788 : if(cplex_fac==2)then
789 : do ii=1,ndgxdtfac
790 : dgxdtfac_(2,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
791 : dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor)
792 : end do
793 : else
794 : do ii=1,ndgxdtfac
795 : ic = cplex_dgxdt(ii) ; jc = 3-ic
796 : dgxdtfac_(ic,ii,ilmn)=zero
797 : if(ic==1)then
798 : dgxdtfac_(jc,ii,ilmn)=-scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
799 : else
800 : dgxdtfac_(jc,ii,ilmn)= scale*dgxdtfac(1,ii,ilmn,ia,ispinor)
801 : end if
802 : end do
803 : end if
804 : end if
805 : end do
806 : !$OMP END DO
807 : end if
808 : if (choice==54.or.choice==8.or.choice==81.or.choice==33) then
809 : !$OMP DO
810 : do ilmn=1,nlmn
811 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
812 : scale=wt;if (il>1) scale=-scale
813 : if (parity) then
814 : if(cplex_fac==2)then
815 : d2gxdtfac_(1:cplex_fac,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac(1:cplex_fac,1:nd2gxdtfac,ilmn,ia,ispinor)
816 : else
817 : do ii=1,nd2gxdtfac
818 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
819 : d2gxdtfac_(ic,ii,ilmn)=scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
820 : d2gxdtfac_(jc,ii,ilmn)=zero
821 : end do
822 : end if
823 : else
824 : if(cplex_fac==2)then
825 : do ii=1,nd2gxdtfac
826 : d2gxdtfac_(1,ii,ilmn)= scale*d2gxdtfac(2,ii,ilmn,ia,ispinor)
827 : d2gxdtfac_(2,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
828 : end do
829 : else
830 : do ii=1,nd2gxdtfac
831 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
832 : d2gxdtfac_(ic,ii,ilmn)=zero
833 : if(ic==1)then
834 : d2gxdtfac_(jc,ii,ilmn)=-scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
835 : else
836 : d2gxdtfac_(jc,ii,ilmn)= scale*d2gxdtfac(1,ii,ilmn,ia,ispinor)
837 : end if
838 : end do
839 : end if
840 : end if
841 : end do
842 : !$OMP END DO
843 : end if
844 : !$OMP END PARALLEL
845 : end if
846 :
847 : ! Scale gxfac_sij with 4pi/sqr(omega).(-i)^l
848 : if (paw_opt>=3) then
849 : !$OMP PARALLEL PRIVATE(ilmn,il,parity,scale,ii,ic,jc)
850 : !$OMP DO
851 : do ilmn=1,nlmn
852 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
853 : scale=wt;if (il>1) scale=-scale
854 : if (parity) then
855 : gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor)
856 : if (cplex==1) gxfacs_(2,ilmn)=zero
857 : else
858 : gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor)
859 : if (cplex==2) then
860 : gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor)
861 : else
862 : gxfacs_(1,ilmn)=zero
863 : end if
864 : end if
865 : end do
866 : !$OMP END DO
867 : if (choice>1) then
868 : !$OMP DO
869 : do ilmn=1,nlmn
870 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
871 : scale=wt;if (il>1) scale=-scale
872 : if (parity) then
873 : if(cplex==2)then
874 : dgxdtfacs_(1:cplex,1:ndgxdtfac,ilmn)=scale*dgxdtfac_sij(1:cplex,1:ndgxdtfac,ilmn,ia,ispinor)
875 : else
876 : do ii=1,ndgxdtfac
877 : ic = cplex_dgxdt(ii) ; jc = 3-ic
878 : dgxdtfacs_(ic,ii,ilmn)=scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
879 : dgxdtfacs_(jc,ii,ilmn)=zero
880 : end do
881 : end if
882 : else
883 : if(cplex==2)then
884 : do ii=1,ndgxdtfac
885 : dgxdtfacs_(2,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
886 : dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor)
887 : end do
888 : else
889 : do ii=1,ndgxdtfac
890 : ic = cplex_dgxdt(ii) ; jc = 3-ic
891 : dgxdtfacs_(ic,ii,ilmn)=zero
892 : if(ic==1)then
893 : dgxdtfacs_(jc,ii,ilmn)=-scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
894 : else
895 : dgxdtfacs_(jc,ii,ilmn)= scale*dgxdtfac_sij(1,ii,ilmn,ia,ispinor)
896 : end if
897 : end do
898 : end if
899 : end if
900 : end do
901 : !$OMP END DO
902 : end if
903 : if (choice==54.or.choice==8.or.choice==81) then
904 : !$OMP DO
905 : do ilmn=1,nlmn
906 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
907 : scale=wt;if (il>1) scale=-scale
908 : if (parity) then
909 : if(cplex==2)then
910 : d2gxdtfacs_(1:cplex,1:nd2gxdtfac,ilmn)=scale*d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,ilmn,ia,ispinor)
911 : else
912 : do ii=1,nd2gxdtfac
913 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
914 : d2gxdtfacs_(ic,ii,ilmn)=scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
915 : d2gxdtfacs_(jc,ii,ilmn)=zero
916 : end do
917 : end if
918 : else
919 : if(cplex==2)then
920 : do ii=1,nd2gxdtfac
921 : d2gxdtfacs_(1,ii,ilmn)= scale*d2gxdtfac_sij(2,ii,ilmn,ia,ispinor)
922 : d2gxdtfacs_(2,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
923 : end do
924 : else
925 : do ii=1,nd2gxdtfac
926 : ic = cplex_d2gxdt(ii) ; jc = 3-ic
927 : d2gxdtfacs_(ic,ii,ilmn)=zero
928 : if(ic==1)then
929 : d2gxdtfacs_(jc,ii,ilmn)=-scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
930 : else
931 : d2gxdtfacs_(jc,ii,ilmn)= scale*d2gxdtfac_sij(1,ii,ilmn,ia,ispinor)
932 : end if
933 : end do
934 : end if
935 : end if
936 : end do
937 : !$OMP END DO
938 : end if
939 : !$OMP END PARALLEL
940 : end if
941 :
942 : ! Compute <g|Vnl|c> (or derivatives) for each plane wave:
943 : if (paw_opt/=3) then
944 : !$OMP PARALLEL PRIVATE(ipw,ilmn,ffnl_dir1)
945 :
946 : ! ------
947 : if (choice==1) then ! <g|Vnl|c>
948 : !$OMP DO
949 : do ipw=1,npw
950 : ztab(ipw)=czero
951 : do ilmn=1,nlmn
952 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
953 : end do
954 : end do
955 : !$OMP END DO
956 :
957 : ! ------
958 : else if (choice==2) then ! derivative w.r.t. atm. pos
959 : !$OMP DO
960 : do ipw=1,npw
961 : ztab(ipw)=czero
962 : do ilmn=1,nlmn
963 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp)
964 : end do
965 : ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw)
966 : do ilmn=1,nlmn
967 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
968 : end do
969 : end do
970 : !$OMP END DO
971 :
972 : ! ------
973 : else if (choice==22) then ! mixed derivative w.r.t. atm. pos and q vector (at q=0)
974 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+qdir
975 : !$OMP DO
976 : do ipw=1,npw
977 : ztab(ipw)=czero
978 : if (idir==qdir) then
979 : do ilmn=1,nlmn
980 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
981 : end do
982 : end if
983 : do ilmn=1,nlmn
984 : ztab(ipw)=ztab(ipw)+kpg(ipw,idir)*ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
985 : ztab(ipw)=ztab(ipw)-ffnl(ipw,ffnl_dir1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
986 : end do
987 : ztab(ipw)=ztab(ipw)*two_pi
988 : end do
989 : !$OMP END DO
990 :
991 : ! ------
992 : else if (choice==25) then ! mixed derivative w.r.t. atm. pos and thwo q vectors (at q=0)
993 : !Use same notation as the notes for clarity
994 : ialpha=nalpha(idir)
995 : idelta=nbeta(idir)
996 : igamma=qdir
997 : idelgam=gamma(idelta,igamma)
998 : !$OMP DO
999 : do ipw=1,npw
1000 : ztab(ipw)=czero
1001 : if (ialpha==igamma) then
1002 : do ilmn=1,nlmn
1003 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1004 : end do
1005 : end if
1006 : if (ialpha==idelta) then
1007 : do ilmn=1,nlmn
1008 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1009 : end do
1010 : end if
1011 : do ilmn=1,nlmn
1012 : ztab(ipw)=ztab(ipw)+kpg(ipw,ialpha)*ffnl(ipw,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1013 : ztab(ipw)=ztab(ipw)-ffnl(ipw,4+idelgam,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
1014 : end do
1015 : ztab(ipw)=ztab(ipw)*two_pi
1016 : end do
1017 : !$OMP END DO
1018 : ! ------
1019 : else if (choice==3) then ! derivative w.r.t. strain
1020 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1021 : if (idir<=3) then
1022 : !$OMP DO
1023 : do ipw=1,npw
1024 : ztab(ipw)=czero
1025 : do ilmn=1,nlmn
1026 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) &
1027 : & *cmplx(dgxdtfac_(1,1,ilmn)-gxfac_(1,ilmn),dgxdtfac_(2,1,ilmn)-gxfac_(2,ilmn),kind=dp) &
1028 : & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1029 : end do
1030 : end do
1031 : !$OMP END DO
1032 : else
1033 : !$OMP DO
1034 : do ipw=1,npw
1035 : ztab(ipw)=czero
1036 : do ilmn=1,nlmn
1037 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) &
1038 : & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1039 : end do
1040 : end do
1041 : !$OMP END DO
1042 : end if
1043 :
1044 : ! ------
1045 : else if (choice==33) then ! mixed derivative w.r.t. strain and q vector (at q=0)
1046 : !Use same notation as the notes for clarity
1047 : ibeta=nalpha(idir)
1048 : idelta=nbeta(idir)
1049 : igamma=qdir
1050 : idelgam=gamma(idelta,igamma)
1051 : !$OMP DO
1052 : do ipw=1,npw
1053 : ztab(ipw)=czero
1054 : if (ibeta==igamma) then
1055 : do ilmn=1,nlmn
1056 : ztab(ipw)=ztab(ipw)+onehalf*ffnl(ipw,1+idelta,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1057 : ztab(ipw)=ztab(ipw)+half*ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp)
1058 : end do
1059 : end if
1060 : if (ibeta==idelta) then
1061 : do ilmn=1,nlmn
1062 : ztab(ipw)=ztab(ipw)+onehalf*ffnl(ipw,1+igamma,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1063 : ztab(ipw)=ztab(ipw)+half*ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
1064 : end do
1065 : end if
1066 : do ilmn=1,nlmn
1067 : ztab(ipw)=ztab(ipw)+kpg(ipw,ibeta)*ffnl(ipw,4+idelgam,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1068 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1+idelta,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)
1069 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1+igamma,ilmn)*cmplx(d2gxdtfac_(1,2,ilmn),d2gxdtfac_(2,2,ilmn),kind=dp)
1070 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfac_(1,3,ilmn),d2gxdtfac_(2,3,ilmn),kind=dp)
1071 : end do
1072 : ztab(ipw)=ztab(ipw)*two_piinv
1073 : end do
1074 : !$OMP END DO
1075 :
1076 : ! ------
1077 : else if (choice==5) then ! full derivative w.r.t. k
1078 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1079 : !$OMP DO
1080 : do ipw=1,npw
1081 : ztab(ipw)=czero
1082 : do ilmn=1,nlmn
1083 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp) &
1084 : & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1085 : end do
1086 : end do
1087 : !$OMP END DO
1088 :
1089 : ! ------
1090 : else if (choice==51) then ! right derivative: <G|p>V<dp/dk|psi>
1091 : !$OMP DO
1092 : do ipw=1,npw
1093 : ztab(ipw)=czero
1094 : do ilmn=1,nlmn
1095 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)
1096 : end do
1097 : end do
1098 : !$OMP END DO
1099 :
1100 : ! ------
1101 : else if (choice==52) then ! left derivative: <G|dp/dk>V<p|psi>
1102 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1103 : !$OMP DO
1104 : do ipw=1,npw
1105 : ztab(ipw)=czero
1106 : do ilmn=1,nlmn
1107 : ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfac_(1,ilmn),gxfac_(2,ilmn),kind=dp)
1108 : end do
1109 : end do
1110 : !$OMP END DO
1111 :
1112 : ! ------
1113 : else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>V<dp/dk_(idir+2)|psi>
1114 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=(mod(idir,3)+1)+1
1115 : !$OMP DO
1116 : do ipw=1,npw
1117 : ztab(ipw)=czero
1118 : do ilmn=1,nlmn
1119 : il=indlmn(1,ilmn)
1120 : ztab(ipw)=ztab(ipw) &
1121 : & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(dgxdtfac_(1,2,ilmn),dgxdtfac_(2,2,ilmn),kind=dp)*CONJG(j_dpc**il)*four_pi
1122 : end do
1123 : end do
1124 : !$OMP END DO
1125 :
1126 : ! ------
1127 : else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k
1128 : !$OMP DO
1129 : do ipw=1,npw
1130 : ztab(ipw)=czero
1131 : do ilmn=1,nlmn
1132 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfac_(2,1,ilmn),-dgxdtfac_(1,1,ilmn),kind=dp)
1133 : end do
1134 : ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw)
1135 : do ilmn=1,nlmn
1136 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)
1137 : end do
1138 : end do
1139 : !$OMP END DO
1140 :
1141 : ! ------
1142 : else if (choice==8) then ! full second order derivative w.r.t. k
1143 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1144 : ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
1145 : ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
1146 : !$OMP DO
1147 : do ipw=1,npw
1148 : ztab(ipw)=czero
1149 : do ilmn=1,nlmn
1150 : ztab(ipw)=ztab(ipw) &
1151 : & +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2, ilmn),kind=dp)&
1152 : & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
1153 : & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,2,ilmn), dgxdtfac_(2,2,ilmn),kind=dp)&
1154 : & +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
1155 : end do
1156 : end do
1157 : !$OMP END DO
1158 :
1159 : ! ------
1160 : else if (choice==81) then
1161 : ! partial second order derivative w.r.t. k
1162 : ! full derivative w.r.t. k1, right derivative w.r.t. k2
1163 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1164 : ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
1165 : !$OMP DO
1166 : do ipw=1,npw
1167 : ztab(ipw)=czero
1168 : do ilmn=1,nlmn
1169 : ztab(ipw)=ztab(ipw) &
1170 : & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfac_(1,1,ilmn),d2gxdtfac_(2,1,ilmn),kind=dp)&
1171 : & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfac_(1,1,ilmn), dgxdtfac_(2,1,ilmn),kind=dp)
1172 : end do
1173 : end do
1174 : !$OMP END DO
1175 :
1176 : ! ------
1177 : else
1178 : !$OMP WORKSHARE
1179 : ztab(:)=czero
1180 : !$OMP END WORKSHARE
1181 : end if
1182 :
1183 :
1184 : ! ------
1185 : !$OMP DO
1186 : do ipw=1,npw
1187 : ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp)
1188 : vect(1,ipw+ipwshft)=vect(1,ipw+ipwshft)+real(ztab(ipw))
1189 : vect(2,ipw+ipwshft)=vect(2,ipw+ipwshft)+aimag(ztab(ipw))
1190 : end do
1191 : !$OMP END DO
1192 :
1193 : !$OMP END PARALLEL
1194 : end if
1195 :
1196 : ! Compute <g|S|c> (or derivatives) for each plane wave:
1197 : if (paw_opt>=3) then
1198 : !$OMP PARALLEL PRIVATE(ilmn,ipw,ffnl_dir1)
1199 :
1200 : ! ------
1201 : if (choice==1) then ! <g|S|c>
1202 : !$OMP DO
1203 : do ipw=1,npw
1204 : ztab(ipw)=czero
1205 : do ilmn=1,nlmn
1206 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1207 : end do
1208 : end do
1209 : !$OMP END DO
1210 :
1211 : ! ------
1212 : else if (choice==2) then ! derivative w.r.t. atm. pos
1213 : !$OMP DO
1214 : do ipw=1,npw
1215 : ztab(ipw)=czero
1216 : do ilmn=1,nlmn
1217 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp)
1218 : end do
1219 : ztab(ipw)=two_pi*kpg(ipw,idir)*ztab(ipw)
1220 : do ilmn=1,nlmn
1221 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
1222 : end do
1223 : end do
1224 : !$OMP END DO
1225 :
1226 : ! ------
1227 : else if (choice==3) then ! derivative w.r.t. strain
1228 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1229 : if (idir<=3) then
1230 : !$OMP DO
1231 : do ipw=1,npw
1232 : ztab(ipw)=czero
1233 : do ilmn=1,nlmn
1234 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn) &
1235 : & *cmplx(dgxdtfacs_(1,1,ilmn)-gxfacs_(1,ilmn),dgxdtfacs_(2,1,ilmn)-gxfacs_(2,ilmn),kind=dp)&
1236 : & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1237 : end do
1238 : end do
1239 : !$OMP END DO
1240 : else
1241 : !$OMP DO
1242 : do ipw=1,npw
1243 : ztab(ipw)=czero
1244 : do ilmn=1,nlmn
1245 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) &
1246 : & -ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1247 : end do
1248 : end do
1249 : !$OMP END DO
1250 : end if
1251 :
1252 : ! ------
1253 : else if (choice==5) then ! full derivative w.r.t. k
1254 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1255 : !$OMP DO
1256 : do ipw=1,npw
1257 : ztab(ipw)=czero
1258 : do ilmn=1,nlmn
1259 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp) &
1260 : & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1261 : end do
1262 : end do
1263 : !$OMP END DO
1264 :
1265 : ! ------
1266 : else if (choice==51) then ! right derivative: <G|p>S<dp/dk|psi>
1267 : !$OMP DO
1268 : do ipw=1,npw
1269 : ztab(ipw)=czero
1270 : do ilmn=1,nlmn
1271 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)
1272 : end do
1273 : end do
1274 : !$OMP END DO
1275 :
1276 : ! ------
1277 : else if (choice==52) then ! left derivative: <G|dp/dk>S<p|psi>
1278 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
1279 : !$OMP DO
1280 : do ipw=1,npw
1281 : ztab(ipw)=czero
1282 : do ilmn=1,nlmn
1283 : ztab(ipw)=ztab(ipw)+ffnl(ipw,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
1284 : end do
1285 : end do
1286 : !$OMP END DO
1287 :
1288 : ! ------
1289 : else if (choice==53) then ! twist derivative: <G|dp/dk_(idir+1)>S<dp/dk_(idir+2)|psi>
1290 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=(mod(idir,3)+1)+1
1291 : !$OMP DO
1292 : do ipw=1,npw
1293 : ztab(ipw)=czero
1294 : do ilmn=1,nlmn
1295 : il=indlmn(1,ilmn)
1296 : ztab(ipw)=ztab(ipw) &
1297 : & +ffnl(ipw,ffnl_dir1,ilmn)*cmplx(dgxdtfacs_(1,2,ilmn),dgxdtfacs_(2,2,ilmn),kind=dp)*CONJG(j_dpc**il)*four_pi
1298 : end do
1299 : end do
1300 : !$OMP END DO
1301 :
1302 : ! ------
1303 : else if (choice==54) then ! mixed derivative w.r.t. atm. pos and k
1304 : !$OMP DO
1305 : do ipw=1,npw
1306 : ztab(ipw)=czero
1307 : do ilmn=1,nlmn
1308 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(dgxdtfacs_(2,1,ilmn),-dgxdtfacs_(1,1,ilmn),kind=dp)
1309 : end do
1310 : ztab(ipw)=two_pi*kpg(ipw,idir1(idir))*ztab(ipw)
1311 : do ilmn=1,nlmn
1312 : ztab(ipw)=ztab(ipw)+ffnl(ipw,1,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)
1313 : end do
1314 : end do
1315 : !$OMP END DO
1316 :
1317 : ! ------
1318 : else if (choice==8) then ! full second order derivative w.r.t. k
1319 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1320 : ffnl_dir(1)=idir1(idir); ffnl_dir(2)=idir2(idir)
1321 : ffnl_dir(3) = gamma(ffnl_dir(1),ffnl_dir(2))
1322 : !$OMP DO
1323 : do ipw=1,npw
1324 : ztab(ipw)=czero
1325 : do ilmn=1,nlmn
1326 : ztab(ipw)=ztab(ipw) &
1327 : & +ffnl(ipw,4+ffnl_dir(3),ilmn)*cmplx( gxfacs_(1, ilmn), gxfacs_(2, ilmn),kind=dp)&
1328 : & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
1329 : & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,2,ilmn), dgxdtfacs_(2,2,ilmn),kind=dp)&
1330 : & +ffnl(ipw,1+ffnl_dir(2),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
1331 : end do
1332 : end do
1333 : !$OMP END DO
1334 :
1335 : ! ------
1336 : else if (choice==81) then
1337 : ! partial second order derivative w.r.t. k
1338 : ! full derivative w.r.t. k1, right derivative w.r.t. k2
1339 : !idir= (xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9)
1340 : ffnl_dir(1)=1; if(dimffnl>2) ffnl_dir(1)=idir1(idir)
1341 : !$OMP DO
1342 : do ipw=1,npw
1343 : ztab(ipw)=czero
1344 : do ilmn=1,nlmn
1345 : ztab(ipw)=ztab(ipw) &
1346 : & +ffnl(ipw,1 ,ilmn)*cmplx(d2gxdtfacs_(1,1,ilmn),d2gxdtfacs_(2,1,ilmn),kind=dp)&
1347 : & +ffnl(ipw,1+ffnl_dir(1),ilmn)*cmplx( dgxdtfacs_(1,1,ilmn), dgxdtfacs_(2,1,ilmn),kind=dp)
1348 : end do
1349 : end do
1350 : !$OMP END DO
1351 :
1352 : ! ------
1353 : else
1354 : !$OMP WORKSHARE
1355 : ztab(:)=czero
1356 : !$OMP END WORKSHARE
1357 : end if
1358 :
1359 :
1360 : ! ------
1361 : ! The OMP WORKSHARE directive doesn't have a good performance with Intel Compiler
1362 : ! !$OMP WORKSHARE
1363 : ! ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
1364 : ! !$OMP END WORKSHARE
1365 : ! !$OMP WORKSHARE
1366 : ! svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
1367 : ! svect(2,1+ipwshft:npw+ipwshft)=svect(2,1+ipwshft:npw+ipwshft)+aimag(ztab(:))
1368 : ! !$OMP END WORKSHARE
1369 : !$OMP DO
1370 : do ipw=1,npw
1371 : ztab(ipw)=ztab(ipw)*cmplx(ph3d(1,ipw,iaph3d),-ph3d(2,ipw,iaph3d),kind=dp)
1372 : svect(1,ipw+ipwshft)=svect(1,ipw+ipwshft)+real(ztab(ipw))
1373 : svect(2,ipw+ipwshft)=svect(2,ipw+ipwshft)+aimag(ztab(ipw))
1374 : end do
1375 : !$OMP END DO
1376 : !$OMP END PARALLEL
1377 : end if
1378 :
1379 : ! End loop on atoms
1380 : end do
1381 : ! End loop on spinors
1382 : end do
1383 :
1384 : ! ==========================================================================
1385 : end if
1386 :
1387 35333413 : ABI_FREE(ztab)
1388 :
1389 35333413 : if (paw_opt/=3) then
1390 31993145 : ABI_FREE(gxfac_)
1391 31993145 : if (choice>1) then
1392 3943262 : ABI_FREE(dgxdtfac_)
1393 : end if
1394 31993145 : if (choice==54.or.choice==8.or.choice==81.or.choice==33) then
1395 453189 : ABI_FREE(d2gxdtfac_)
1396 : end if
1397 : end if
1398 35333413 : if (paw_opt>=3) then
1399 18185405 : ABI_FREE(gxfacs_)
1400 18185405 : if (choice>1) then
1401 1700061 : ABI_FREE(dgxdtfacs_)
1402 : end if
1403 18185405 : if (choice==54.or.choice==8.or.choice==81) then
1404 32973 : ABI_FREE(d2gxdtfacs_)
1405 : end if
1406 : end if
1407 :
1408 : DBG_EXIT("COLL")
1409 :
1410 : #if !defined HAVE_OPENMP
1411 : !Fake use of unused variable
1412 : if (.false.) write(std_out,*) ipw
1413 : #endif
1414 :
1415 35333413 : end subroutine opernlb_ylm
1416 : !!***
1417 :
1418 : end module m_opernlb_ylm
1419 : !!***
|
