Branch data 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 : 36345094 : subroutine opernlb_ylm(choice,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,&
153 : 36345094 : & d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnl,ffnl,gxfac,gxfac_sij,&
154 : 36345094 : & ia3,idir,indlmn,kpg,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpg,nlmn,nloalg,npw,&
155 [ - + - + ]: 36345094 : & 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 : 36345094 : real(dp),allocatable :: d2gxdtfac_(:,:,:),d2gxdtfacs_(:,:,:),dgxdtfac_(:,:,:),dgxdtfacs_(:,:,:),gxfac_(:,:),gxfacs_(:,:)
185 : : ! real(dp),allocatable :: kpg(:,:)
186 : 36345094 : complex(dp),allocatable :: ztab(:)
187 : : ! *************************************************************************
188 : :
189 : : DBG_ENTER("COLL")
190 : :
191 : : !Nothing to do when choice=4, 6 or 23
192 [ + - ]: 36345094 : if (choice==4.or.choice==6.or.choice==23) return
193 : :
194 : : !DDK not compatible with istwkf > 1
195 [ + + + + : 78419244 : 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 : 36345094 : wt=four_pi/sqrt(ucvol)
201 : 36345094 : nthreads=1
202 : : #if defined HAVE_OPENMP
203 : : nthreads=OMP_GET_NUM_THREADS()
204 : : #endif
205 : :
206 [ + + ]: 36345094 : if (paw_opt/=3) then
207 [ + + + - : 98018442 : ABI_MALLOC(gxfac_,(2,nlmn))
+ - + + -
+ - + ]
208 [ + + + + ]: 1060310031 : gxfac_(:,:)=zero
209 [ + + ]: 32672814 : if (choice>1) then
210 [ - + + - : 15884072 : ABI_MALLOC(dgxdtfac_,(2,ndgxdtfac,nlmn))
+ + + - +
- + + - +
- + ]
211 [ + - + + : 162188942 : if(ndgxdtfac>0) dgxdtfac_(:,:,:)=zero
+ + + + ]
212 : : end if
213 [ + + ]: 32672814 : 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 [ + + ]: 36345094 : if (paw_opt>=3) then
219 [ - + + - : 57273258 : ABI_MALLOC(gxfacs_,(2,nlmn))
+ - + - -
+ - + ]
220 [ + + + + ]: 643167091 : gxfacs_(:,:)=zero
221 [ + + ]: 19091086 : if (choice>1) then
222 [ - + + - : 6892980 : ABI_MALLOC(dgxdtfacs_,(2,ndgxdtfac,nlmn))
- + + - +
- + - - +
- + ]
223 [ + - + + : 62596332 : if (ndgxdtfac>0) dgxdtfacs_(:,:,:)=zero
+ + + + ]
224 : : end if
225 [ + + + + ]: 19091086 : 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 [ - + ]: 36345094 : 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 [ - + + - : 109035282 : ABI_MALLOC(ztab,(npw))
+ - + - -
+ - + ]
239 : :
240 : : !==========================================================================
241 : : !========== STANDARD VERSION ==============================================
242 : : !==========================================================================
243 : : if (nthreads==1) then
244 : :
245 : : ! Loop on spinorial components
246 [ + + ]: 74702196 : do ispinor=1,nspinor
247 : 38357102 : ipwshft=(ispinor-1)*npw
248 : :
249 : : ! Loop on atoms (blocking)
250 [ + + ]: 125945859 : do ia=1,nincat
251 [ + + ]: 51243663 : iaph3d=ia;if (nloalg(2)>0) iaph3d=ia+ia3-1
252 : : ! Scale gxfac with 4pi/sqr(omega).(-i)^l
253 [ + + ]: 51243663 : if (paw_opt/=3) then
254 [ + + ]: 532042375 : do ilmn=1,nlmn
255 : 485906324 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
256 [ + + ]: 485906324 : scale=wt;if (il>1) scale=-scale
257 [ + + ]: 532042375 : if (parity) then
258 [ + + ]: 611307135 : gxfac_(1:cplex_fac,ilmn)=scale*gxfac(1:cplex_fac,ilmn,ia,ispinor)
259 [ + + ]: 207071377 : if (cplex_fac==1) gxfac_(2,ilmn)=zero
260 : : else
261 : 278834947 : gxfac_(2,ilmn)=-scale*gxfac(1,ilmn,ia,ispinor)
262 [ + + ]: 278834947 : if (cplex_fac==2) then
263 : 255224651 : gxfac_(1,ilmn)=scale*gxfac(2,ilmn,ia,ispinor)
264 : : else
265 : 23610296 : gxfac_(1,ilmn)=zero
266 : : end if
267 : : end if! parity
268 : : end do ! ilmn
269 [ + + ]: 46136051 : if (choice>1) then
270 [ + + ]: 51980764 : do ilmn=1,nlmn
271 : 46872291 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
272 [ + + ]: 46872291 : scale=wt;if (il>1) scale=-scale
273 [ + + ]: 51980764 : if (parity) then
274 [ + - ]: 20286669 : if(cplex_fac==2)then
275 [ + + + + ]: 93034479 : 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 [ + - ]: 26585622 : if(cplex_fac==2)then
285 [ + + ]: 56644578 : do ii=1,ndgxdtfac
286 : 30058956 : dgxdtfac_(1,ii,ilmn)= scale*dgxdtfac(2,ii,ilmn,ia,ispinor)
287 : 56644578 : 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 [ + + ]: 46136051 : 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 [ + + ]: 51243663 : if (paw_opt>=3) then
341 : :
342 [ + + ]: 354698603 : do ilmn=1,nlmn
343 : 324702594 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
344 [ + + ]: 324702594 : scale=wt;if (il>1) scale=-scale
345 [ + + ]: 354698603 : if (parity) then
346 [ + + ]: 414231912 : gxfacs_(1:cplex,ilmn)=scale*gxfac_sij(1:cplex,ilmn,ia,ispinor)
347 [ + + ]: 141357504 : if (cplex==1) gxfacs_(2,ilmn)=zero
348 : : else
349 : 183345090 : gxfacs_(2,ilmn)=-scale*gxfac_sij(1,ilmn,ia,ispinor)
350 [ + + ]: 183345090 : if (cplex==2) then
351 : 159199238 : gxfacs_(1,ilmn)=scale*gxfac_sij(2,ilmn,ia,ispinor)
352 : : else
353 : 24145852 : gxfacs_(1,ilmn)=zero
354 : : end if
355 : : end if
356 : : end do
357 [ + + ]: 29996009 : if (choice>1) then
358 [ + + ]: 19252100 : do ilmn=1,nlmn
359 : 17279810 : il=mod(indlmn(1,ilmn),4);parity=(mod(il,2)==0)
360 [ + + ]: 17279810 : scale=wt;if (il>1) scale=-scale
361 [ + + ]: 19252100 : if (parity) then
362 [ + - ]: 5385698 : if(cplex==2)then
363 [ + + + + ]: 21591158 : 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 [ + - ]: 11894112 : if(cplex==2)then
373 [ + + ]: 23818500 : do ii=1,ndgxdtfac
374 : 11924388 : dgxdtfacs_(1,ii,ilmn)= scale*dgxdtfac_sij(2,ii,ilmn,ia,ispinor)
375 : 23818500 : 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 [ + + + + ]: 29996009 : 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 [ + + ]: 30272807 : 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 [ + + ]: 51243663 : if (paw_opt/=3) then
430 : :
431 [ + + ]: 14316395262 : ztab(:)=czero
432 : :
433 : : ! ------
434 [ + + ]: 46136051 : if (choice==1) then ! <g|Vnl|c>
435 [ + + ]: 480061611 : do ilmn=1,nlmn
436 [ + + ]: >14066*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 [ + + ]: 46136051 : if (choice==2) then ! derivative w.r.t. atm. pos
442 [ + + ]: 15170090 : do ilmn=1,nlmn
443 [ + + ]: 2081417571 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfac_(2,ilmn),-gxfac_(1,ilmn),kind=dp)
444 : : end do
445 [ + + ]: 223431097 : ztab(:)=two_pi*kpg(:,idir)*ztab(:)
446 [ + + ]: 15170090 : do ilmn=1,nlmn
447 [ + + ]: 2081417571 : 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 [ + + ]: 46136051 : 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 [ + + ]: 46136051 : 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 [ + + ]: 46136051 : 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 [ + + ]: 46136051 : 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 [ + + ]: 46136051 : if (choice==5) then ! full derivative w.r.t. k
537 [ + + ]: 1701532 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
538 [ + + ]: 17643774 : do ilmn=1,nlmn
539 : : ztab(:)=ztab(:)+ffnl(:,1 ,ilmn)*cmplx(dgxdtfac_(1,1,ilmn),dgxdtfac_(2,1,ilmn),kind=dp)&
540 [ + + ]: 2569116768 : & +ffnl(:,ffnl_dir1,ilmn)*cmplx( gxfac_(1, ilmn), gxfac_(2 ,ilmn),kind=dp)
541 : : end do
542 : : end if
543 : :
544 : : ! ------
545 [ - + ]: 46136051 : 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 [ - + ]: 46136051 : 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 [ - + ]: 46136051 : 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 [ - + ]: 46136051 : 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 [ + + ]: 46136051 : 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 [ - + ]: 46136051 : 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 [ + + ]: 14316395262 : ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
609 : :
610 [ + + ]: 14316395262 : vect(1,1+ipwshft:npw+ipwshft)=vect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
611 [ + + ]: 14316395262 : 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 [ + + ]: 89600765 : if (paw_opt>=3) then
618 : :
619 [ + + ]: 11830817174 : ztab(:)=czero
620 : :
621 : : ! ------
622 [ + + ]: 29996009 : if (choice==1) then ! <g|S|c>
623 [ + + ]: 335446503 : do ilmn=1,nlmn
624 [ + + ]: >12028*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 [ + + ]: 29996009 : if (choice==2) then ! derivative w.r.t. atm. pos
630 [ + + ]: 7438624 : do ilmn=1,nlmn
631 [ + + ]: 944639620 : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(gxfacs_(2,ilmn),-gxfacs_(1,ilmn),kind=dp)
632 : : end do
633 [ + + ]: 106451388 : ztab(:)=two_pi*kpg(:,idir)*ztab(:)
634 [ + + ]: 7438624 : do ilmn=1,nlmn
635 [ + + ]: 944639620 : 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 [ + + ]: 29996009 : 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 [ + + ]: 29996009 : if (choice==5) then ! full derivative w.r.t. k
658 [ + + ]: 613976 : ffnl_dir1=2; if(dimffnl>2) ffnl_dir1=1+idir
659 [ + + ]: 6388504 : do ilmn=1,nlmn
660 : : ztab(:)=ztab(:)+ffnl(:,1,ilmn)*cmplx(dgxdtfacs_(1,1,ilmn),dgxdtfacs_(2,1,ilmn),kind=dp)&
661 [ + + ]: 978012744 : & +ffnl(:,ffnl_dir1,ilmn)*cmplx(gxfacs_(1,ilmn),gxfacs_(2,ilmn),kind=dp)
662 : : end do
663 : : end if
664 : :
665 : : ! ------
666 [ + + ]: 29996009 : 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 [ - + ]: 29996009 : 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 [ - + ]: 29996009 : 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 [ + + ]: 29996009 : 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 [ + + ]: 29996009 : 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 [ + + ]: 29996009 : 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 [ + + ]: 11830817174 : ztab(:)=ztab(:)*cmplx(ph3d(1,:,iaph3d),-ph3d(2,:,iaph3d),kind=dp)
730 [ + + ]: 11830817174 : svect(1,1+ipwshft:npw+ipwshft)=svect(1,1+ipwshft:npw+ipwshft)+real(ztab(:))
731 [ + + ]: 11830817174 : 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 [ - + ]: 36345094 : ABI_FREE(ztab)
1388 : :
1389 [ + + ]: 36345094 : if (paw_opt/=3) then
1390 [ - + ]: 32672814 : ABI_FREE(gxfac_)
1391 [ + + ]: 32672814 : if (choice>1) then
1392 [ - + ]: 3971018 : ABI_FREE(dgxdtfac_)
1393 : : end if
1394 [ + + ]: 32672814 : 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 [ + + ]: 36345094 : if (paw_opt>=3) then
1399 [ - + ]: 19091086 : ABI_FREE(gxfacs_)
1400 [ + + ]: 19091086 : if (choice>1) then
1401 [ - + ]: 1723245 : ABI_FREE(dgxdtfacs_)
1402 : : end if
1403 [ + + + + ]: 19091086 : 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 [ - + - + : 36345094 : end subroutine opernlb_ylm
- + - + -
+ - + ]
1416 : : !!***
1417 : :
1418 : : end module m_opernlb_ylm
1419 : : !!***
|
