Branch data Line data Source code
1 : : !!****m* ABINIT/m_opernl
2 : : !! NAME
3 : : !! m_opernl
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, DRH)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_opernl
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : :
28 : : use m_mkffkg, only : mkffkg, dfpt_mkffkg
29 : :
30 : : implicit none
31 : :
32 : : private
33 : : !!***
34 : :
35 : : public :: opernl2
36 : : public :: opernl3
37 : : public :: opernl4a
38 : : public :: opernl4b
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* ABINIT/opernl2
45 : : !! NAME
46 : : !! opernl2
47 : : !!
48 : : !! FUNCTION
49 : : !! Operate with the non-local part of the hamiltonian,
50 : : !! either from reciprocal space to projected quantities (sign=1),
51 : : !! or from projected quantities to reciprocal space (sign=-1)
52 : : !!
53 : : !! INPUTS
54 : : !! if(sign==-1 .and. (choice==2 .or choice==4 .or. choice==5))
55 : : !! dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
56 : : !! or with respect to ddk
57 : : !! if(sign==-1 .and. choice==3)
58 : : !! dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
59 : : !! ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
60 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
61 : : !! ia3=gives the number of the first atom in the subset presently treated
62 : : !! idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
63 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
64 : : !! ispinor=1 or 2, gives the spinorial component of ffnl to be used
65 : : !! istwf_k=option parameter that describes the storage of wfs
66 : : !! itypat = type of atom, needed for ffnl
67 : : !! jproj(nlang)=number of projectors for each angular momentum
68 : : !! kg_k(3,npw)=integer coords of planewaves in basis sphere
69 : : !! kpg_k(npw,npkg)= (k+G) components and related data
70 : : !! kpt(3)=real components of k point in terms of recip. translations
71 : : !! lmnmax=max. number of (l,n) components over all type of psps
72 : : !! matblk=dimension of the array ph3d
73 : : !! mincat= maximum increment of atoms
74 : : !! mlang1 = dimensions for dgxdis1
75 : : !! mlang3 = one of the dimensions of the array gxa
76 : : !! mlang4 = dimension for dgxds
77 : : !! mlang5 = dimensions for dgxdis2
78 : : !! mlang6 = dimension for d2gxds2
79 : : !! mproj=maximum dimension for number of projection operators for each
80 : : !! angular momentum for nonlocal pseudopotential
81 : : !! ndgxdt=second dimension of dgxdt
82 : : !! nffnl=second dimension of ffnl
83 : : !! nincat = number of atoms in the subset here treated
84 : : !! nkpg=second size of array kpg_k
85 : : !! nlang = number of angular momenta to be treated = 1 + highest ang. mom.
86 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
87 : : !! npw = number of plane waves in reciprocal space
88 : : !! ntypat = number of type of atoms, dimension needed for ffnl
89 : : !! sign : if 1, go from reciprocal space to projected scalars,
90 : : !! if -1, go from projected scalars to reciprocal space.
91 : : !! if(sign==1), vect(2*npw)=starting vector in reciprocal space
92 : : !! if(sign==-1) gxa(2,mlang3,nincat,mproj)=modified projected scalars;
93 : : !! NOTE that metric contractions have already been performed on the
94 : : !! arrays gxa if sign=-1
95 : : !! ph3d(2,npw,matblk)=three-dimensional phase factors
96 : : !!
97 : : !! OUTPUT
98 : : !! if(sign==1)
99 : : !! gxa(2,mlang3,mincat,mproj)= projected scalars
100 : : !! if(sign==1 .and. (choice==2 .or choice==4 .or. choice==5 .or. choice==23))
101 : : !! dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
102 : : !! or with respect to ddk
103 : : !! if(sign==1 .and. (choice==3 .or. choice==23))
104 : : !! dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
105 : : !! if(sign==1 .and. choice==6)
106 : : !! dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars
107 : : !! wrt coord. indexed for internal strain
108 : : !! d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars
109 : : !! wrt strain and coord
110 : : !! d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars
111 : : !! wrt strains
112 : : !! if(sign==-1)
113 : : !! vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.
114 : : !!
115 : : !! NOTES
116 : : !! Operate with the non-local part of the hamiltonian for one type of
117 : : !! atom, and within this given type of atom, for a subset
118 : : !! of at most nincat atoms.
119 : : !!
120 : : !! This routine basically replaces getgla (gxa here is the former gla),
121 : : !! except for the calculation of <G|dVnl/dk|C> or strain gradients.
122 : : !!
123 : : !! SOURCE
124 : :
125 : 20868 : subroutine opernl2(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
126 : 20868 : & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
127 : 20868 : & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
128 : : & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
129 [ - + ]: 20868 : & ntypat,ph3d,sign,vect)
130 : :
131 : :
132 : : !Arguments ------------------------------------
133 : : !scalars
134 : : integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk
135 : : integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt
136 : : integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign
137 : : !arrays
138 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
139 : : integer,intent(in) :: nloalg(3)
140 : : real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
141 : : real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk)
142 : : real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj)
143 : : real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
144 : : real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:)
145 : : real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj)
146 : : real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj)
147 : : real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj)
148 : :
149 : : !Local variables-------------------------------
150 : : !scalars
151 : : integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3
152 : : integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw
153 : : integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj,ntens
154 : : real(dp),parameter :: two_pi2=two_pi**2
155 : : character(len=500) :: message
156 : : !arrays
157 : 20868 : integer,allocatable :: parity(:)
158 : : !real(dp) :: tsec(2)
159 : 20868 : real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
160 : :
161 : : ! *************************************************************************
162 : :
163 : : !call wrtout(std_out,"in opernl2","COLL")
164 : :
165 : : !mblkpw sets the size of blocks of planewaves
166 : 20868 : mblkpw=NLO_MBLKPW
167 : :
168 : : !Get the actual maximum number of projectors
169 [ + + ]: 104340 : mmproj=maxval(indlmn(3,:,itypat))
170 : :
171 : : !Initialisation before blocking on the plane waves
172 : :
173 [ + + ]: 20868 : if (sign==1) then
174 : : ! Put projected scalars to zero
175 [ + + + + : 833214 : gxa(:,:,:,1:mmproj)=0.0d0
+ + + + ]
176 [ + + + + : 39386 : if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
+ + + + +
+ + + ]
177 [ + + + + : 213978 : if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
+ + + + +
+ ]
178 [ + + ]: 10474 : if (choice==6) then
179 [ + + + + : 3360 : dgxdis(:,:,:,1:mmproj)=0.0d0
+ + + + ]
180 [ + + + + : 7200 : d2gxdis(:,:,:,1:mmproj)=0.0d0
+ + + + ]
181 [ + + + + : 9696 : d2gxds2(:,:,:,1:mmproj)=0.0d0
+ + + + ]
182 : : end if
183 : : end if
184 : :
185 : : !Set up dimension of kpgx and allocate
186 : : !ntens sets the maximum number of independent tensor components
187 : : !over all allowed angular momenta; need 20 for spdf for tensors
188 : : !up to rank 3; to handle stress tensor, need up to rank 5
189 : 20868 : ntens=1
190 [ - + - - ]: 20868 : if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4
191 [ + + + + ]: 20868 : if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
192 [ + + + + : 20868 : if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
+ - ]
193 [ + + + + : 20868 : if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
+ + ]
194 [ + + + + : 20868 : if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
+ + + - ]
195 [ + + + + ]: 20868 : if(choice==6 .and. nlang>=3)ntens=84
196 [ + + + + ]: 20868 : if(choice==6 .and. nlang==4)ntens=120
197 : :
198 : : !Set up second dimension of ffkg array, and allocate
199 : 20868 : nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
200 [ + + ]: 80684 : do ilang=1,nlang
201 : : ! Get the number of projectors for that angular momentum
202 : 59816 : nproj=jproj(ilang)
203 : : ! If there is a non-local part, accumulate the number of vectors needed
204 : : ! The variables ilang below are the number of independent tensors of
205 : : ! various ranks, the variable names being more historical than logical.
206 : : ! ilang2=number of rank ilang-1
207 : : ! ilang3=number of rank ilang+1
208 : : ! ilang4=number of rank ilang
209 : : ! ilang5=number of rank ilang+2
210 : : ! ilang6=number of rank ilang+3
211 [ + + ]: 80684 : if(nproj>0)then
212 : 52170 : ilang2=(ilang*(ilang+1))/2
213 : 52170 : nffkge=nffkge+nproj*ilang2
214 [ - + ]: 52170 : if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
215 [ + + ]: 52170 : if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2
216 [ + + ]: 52170 : if(choice==3 .or. choice==6 .or. choice==23)then
217 : 6120 : ilang3=((ilang+2)*(ilang+3))/2
218 : 6120 : nffkgs=nffkgs+nproj*ilang3
219 : : end if
220 [ + + ]: 52170 : if(choice==6)then
221 : 80 : ilang4=((ilang+1)*(ilang+2))/2
222 : 80 : ilang5=((ilang+3)*(ilang+4))/2
223 : 80 : ilang6=((ilang+4)*(ilang+5))/2
224 : 80 : nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6)
225 : : end if
226 : : end if
227 : : end do
228 : 20868 : nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk
229 : :
230 : : !Loop on subsets of plane waves (blocking)
231 : : !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a sigfault with nthreads==1
232 : : !Feb 16 2012: The code does not crash anymore but it's not efficient.
233 : : !
234 : : !!$OMP PARALLEL DEFAULT(PRIVATE) &
235 : : !!$OMP SHARED(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) &
236 : : !!$OMP SHARED(ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) &
237 : : !!$OMP SHARED(jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) &
238 : : !!$OMP SHARED(mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) &
239 : : !!$OMP SHARED(ntypat,ph3d,sign,vect) &
240 : : !!$OMP SHARED(mblkpw,nffkg,nffkgd,nffkge,nffkgs,ntens,mmproj)
241 : :
242 [ - + + - : 62604 : ABI_MALLOC(ffkg,(mblkpw,nffkg))
+ - + - -
+ - + ]
243 [ - + + - : 62604 : ABI_MALLOC(parity,(nffkg))
+ - + - -
+ - + ]
244 [ - + + - : 62604 : ABI_MALLOC(kpgx,(mblkpw,ntens))
+ - + - -
+ - + ]
245 [ - + + - : 62604 : ABI_MALLOC(scalars,(2,nffkg))
+ - + - -
+ - + ]
246 [ - + ]: 20868 : ABI_MALLOC(teffv,(2,mblkpw))
247 : : !!$OMP DO
248 [ + - + + ]: 112272 : do ipw1=1,npw,mblkpw
249 : :
250 : 91404 : ipw2=min(npw,ipw1+mblkpw-1)
251 : 91404 : nincpw=ipw2-ipw1+1
252 : :
253 : : ! call timab(74+choice,1,tsec)
254 : :
255 : : ! Initialize kpgx array related to tensors defined below
256 : : call mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,kg_k,&
257 : : & kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
258 : 91404 : & npw,ntens,ntypat,parity)
259 : :
260 : : ! call timab(74+choice,2,tsec)
261 : :
262 : : ! Now treat the different signs
263 [ + + ]: 112272 : if (sign==1) then
264 : :
265 [ + + ]: 91788 : do ia=1,nincat
266 : :
267 : : ! Compute the shift eventually needed to get the phases in ph3d
268 : 45894 : iaph3d=ia
269 [ + + ]: 45894 : if(nloalg(2)>0)iaph3d=ia+ia3-1
270 : :
271 : : ! ******* Entering the first time-consuming part of the routine *******
272 : :
273 : : ! Multiply by the phase factor
274 : : ! This allows to be left with only real operations,
275 : : ! that are performed in the most inner loops
276 : 45894 : ig=ipw1
277 [ + + ]: 7776088 : do ipw=1,nincpw
278 : 7730194 : teffv(1,ipw)=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
279 : 7730194 : teffv(2,ipw)=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
280 : 7776088 : ig=ig+1
281 : : end do
282 : :
283 [ + + ]: 645993 : do iffkg=1,nffkg
284 : 600099 : scalars(1,iffkg)=0.0d0
285 : 600099 : scalars(2,iffkg)=0.0d0
286 [ + + ]: 98904460 : do ipw=1,nincpw
287 : 98258467 : scalars(1,iffkg)=scalars(1,iffkg)+teffv(1,ipw)*ffkg(ipw,iffkg)
288 : 98858566 : scalars(2,iffkg)=scalars(2,iffkg)+teffv(2,ipw)*ffkg(ipw,iffkg)
289 : : end do
290 : : end do
291 : :
292 : : ! ******* Leaving the critical part *********************************
293 : :
294 : : ! DEBUG
295 : : ! write(std_out,*)' opernl2 : scalars'
296 : : ! do iffkg=1,10
297 : : ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
298 : : ! end do
299 : : ! stop
300 : : ! ENDDEBUG
301 : :
302 [ + + ]: 45894 : if(istwf_k>=2)then
303 : : ! Impose parity of resulting scalar (this operation could be
304 : : ! replaced by direct saving of CPU time in the preceeding section)
305 [ + + ]: 112266 : do iffkg=1,nffkg
306 : 112266 : scalars(parity(iffkg),iffkg)=0.0d0
307 : : end do
308 : : end if
309 : :
310 : 45894 : iffkg=0
311 : 45894 : iffkgs=nffkge+nffkgd
312 : 45894 : iffkgs2=nffkge+nffkgs
313 : 45894 : iffkgk=nffkge*2
314 [ + + ]: 218040 : do ilang=1,nlang
315 : 126252 : nproj=jproj(ilang)
316 [ + + ]: 172146 : if(nproj>0)then
317 : : ! ilang2 is the number of independent tensor components
318 : : ! for symmetric tensor of rank ilang-1
319 : 114735 : ilang2=(ilang*(ilang+1))/2
320 : :
321 : : ! Loop over projectors
322 [ + + ]: 263760 : do iproj=1,nproj
323 : : ! Multiply by the k+G factors (tensors of various rank)
324 [ + + ]: 608052 : do ii=1,ilang2
325 : : ! Get the starting address for the relevant tensor
326 : 459027 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
327 : 459027 : iffkg=iffkg+1
328 : : ! !$OMP CRITICAL (OPERNL2_1)
329 : 459027 : gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
330 : 459027 : gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg)
331 : : ! !$OMP END CRITICAL (OPERNL2_1)
332 : : ! Now, compute gradients, if needed.
333 [ + + + + ]: 459027 : if ((choice==2.or.choice==23) .and. ndgxdt==3) then
334 [ + + ]: 11328 : do mu=1,3
335 : 8496 : jffkg=nffkge+(iffkg-1)*3+mu
336 : : ! Pay attention to the use of reals and imaginary parts here ...
337 : : ! !$OMP CRITICAL (OPERNL2_2)
338 : 8496 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
339 : 11328 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
340 : : ! !$OMP END CRITICAL (OPERNL2_2)
341 : : end do
342 : : end if
343 [ + + + - ]: 459027 : if (choice==2 .and. ndgxdt==1) then
344 : 6048 : jffkg=nffkge+iffkg
345 : : ! Pay attention to the use of reals and imaginary parts here ...
346 : : ! !$OMP CRITICAL (OPERNL2_3)
347 : 6048 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg)
348 : 6048 : dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg)
349 : : ! !$OMP END CRITICAL (OPERNL2_3)
350 : : end if
351 [ - + ]: 608052 : if (choice==4) then
352 [ # # ]: 0 : do mu=1,3
353 : 0 : jffkg=nffkge+(iffkg-1)*9+mu
354 : : ! Pay attention to the use of reals and imaginary parts here ...
355 : : ! !$OMP CRITICAL (OPERNL2_4)
356 : 0 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
357 : 0 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
358 : : ! !$OMP END CRITICAL (OPERNL2_4)
359 : : end do
360 [ # # ]: 0 : do mu=4,9
361 : 0 : jffkg=nffkge+(iffkg-1)*9+mu
362 : : ! Pay attention to the use of reals and imaginary parts here ...
363 : : ! Also, note the multiplication by (2 pi)**2
364 : : ! !$OMP CRITICAL (OPERNL2_5)
365 : 0 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg)
366 : 0 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg)
367 : : ! !$OMP END CRITICAL (OPERNL2_5)
368 : : end do
369 : : end if
370 : : ! End loop on ii=1,ilang2
371 : : end do
372 : :
373 [ + + ]: 149025 : if (choice==3 .or. choice==6 .or. choice==23) then
374 : : ! Compute additional tensors related to strain gradients
375 : : ! ilang3 is number of unique tensor components of rank ilang+1
376 : 10128 : ilang3=((ilang+2)*(ilang+3))/2
377 : 10128 : jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
378 : : ! Compute strain gradient tensor components
379 [ + + ]: 125328 : do ii=1,ilang3
380 : : ! Note that iffkgs is also used by ddk and 2nd derivative parts
381 : 115200 : iffkgs=iffkgs+1
382 : 115200 : jj=ii+jjs
383 : : ! !$OMP CRITICAL (OPERNL2_6)
384 : 115200 : dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
385 : 125328 : dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs)
386 : : ! !$OMP END CRITICAL (OPERNL2_6)
387 : : end do
388 : : end if
389 : :
390 [ + + ]: 149025 : if (choice==6) then
391 : : ! Compute additional tensors related to strain 2nd derivatives
392 : : ! and internal strain derivatives
393 : : ! ilang6 is number of unique tensor components of rank ilang+3
394 : 240 : ilang6=((ilang+4)*(ilang+5))/2
395 : 240 : jjs=((ilang+3)*(ilang+4)*(ilang+5))/6
396 : : ! Compute strain gradient tensor components
397 [ + + ]: 5760 : do ii=1,ilang6
398 : 5520 : iffkgs2=iffkgs2+1
399 : 5520 : jj=ii+jjs
400 : : ! !$OMP CRITICAL (OPERNL2_6)
401 : 5520 : d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2)
402 : 5760 : d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2)
403 : : ! !$OMP END CRITICAL (OPERNL2_6)
404 : : end do
405 : :
406 : : ! ilang4 is number of unique tensor components of rank ilang
407 : 240 : ilang4=((ilang+1)*(ilang+2))/2
408 : 240 : jjs=((ilang)*(ilang+1)*(ilang+2))/6
409 : : ! Compute internal strain gradient tensor components
410 [ + + ]: 2016 : do ii=1,ilang4
411 : 1776 : iffkgs2=iffkgs2+1
412 : 1776 : jj=ii+jjs
413 : : ! !$OMP CRITICAL (OPERNL2_6)
414 : : ! Pay attention to the use of reals and imaginary parts here ...
415 : 1776 : dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2)
416 : 2016 : dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2)
417 : : ! !$OMP END CRITICAL (OPERNL2_6)
418 : : end do
419 : :
420 : : ! ilang5 is number of unique tensor components of rank ilang+2
421 : 240 : ilang5=((ilang+3)*(ilang+4))/2
422 : 240 : jjs=((ilang+2)*(ilang+3)*(ilang+4))/6
423 : : ! Compute internal strain gradient tensor components
424 [ + + ]: 4272 : do ii=1,ilang5
425 : 4032 : iffkgs2=iffkgs2+1
426 : 4032 : jj=ii+jjs
427 : : ! !$OMP CRITICAL (OPERNL2_6)
428 : : ! Pay attention to the use of reals and imaginary parts here ...
429 : 4032 : d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2)
430 : 4272 : d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2)
431 : : ! !$OMP END CRITICAL (OPERNL2_6)
432 : : end do
433 : : end if ! choice==6
434 : :
435 [ - + ]: 263760 : if (choice==5) then
436 : : ! Compute additional tensors related to ddk with ffnl(:,2,...)
437 : 0 : ilangx=(ilang*(ilang+1))/2
438 : 0 : jjs=((ilang-1)*ilang*(ilang+1))/6
439 [ # # ]: 0 : do ii=1,ilangx
440 : : ! Note that iffkgs is also used by strain part
441 : 0 : iffkgs=iffkgs+1
442 : 0 : jj=ii+jjs
443 : : ! !$OMP CRITICAL (OPERNL2_7)
444 : 0 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
445 : 0 : dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs)
446 : : ! !$OMP END CRITICAL (OPERNL2_7)
447 : : end do
448 : : ! Compute additional tensors related to ddk with ffnl(:,1,...)
449 [ # # ]: 0 : if(ilang>=2)then
450 : 0 : ilangx=((ilang-1)*ilang)/2
451 : 0 : jjs=((ilang-2)*(ilang-1)*ilang)/6
452 [ # # ]: 0 : do ii=1,ilangx
453 : 0 : iffkgk=iffkgk+1
454 : 0 : jj=ii+jjs
455 : : ! !$OMP CRITICAL (OPERNL2_8)
456 : 0 : dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
457 : 0 : dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk)
458 : : ! !$OMP END CRITICAL (OPERNL2_8)
459 : : end do
460 : : end if
461 : : end if
462 : :
463 : : ! End projector loop
464 : : end do
465 : :
466 : : ! End condition of non-zero projectors
467 : : end if
468 : :
469 : : ! End angular momentum loop
470 : : end do
471 : :
472 : : ! End loop on atoms
473 : : end do
474 : :
475 [ + - ]: 45510 : else if (sign==-1) then
476 : : ! Application of non-local part from projected scalars
477 : : ! back to reciprocal space ...
478 : : ! [this section merely computes terms which add to <G|Vnl|C>;
479 : : ! nothing here is needed when various gradients are being computed]
480 : :
481 : : ! Loop on atoms
482 [ + + ]: 91020 : do ia=1,nincat
483 : :
484 : : ! Compute the shift eventually needed to get the phases in ph3d
485 : 45510 : iaph3d=ia
486 [ + + ]: 45510 : if(nloalg(2)>0)iaph3d=ia+ia3-1
487 : :
488 : : ! Transfer gxa (and eventually dgxdt) in scalars with different indexing
489 : 45510 : iffkg=0
490 : 45510 : iffkgk=nffkge*2
491 : 45510 : iffkgs=nffkge
492 [ + + ]: 170706 : do ilang=1,nlang
493 : 125196 : nproj=jproj(ilang)
494 [ + + ]: 170706 : if (nproj>0) then
495 : 113775 : ilang2=(ilang*(ilang+1))/2
496 : 113775 : ilang3=((ilang+2)*(ilang+3))/2
497 [ + + ]: 261552 : do iproj=1,nproj
498 [ + + ]: 602964 : do ii=1,ilang2
499 : 455187 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
500 : 455187 : iffkg=iffkg+1
501 [ + + ]: 602964 : if(choice==1 .or. choice==3)then
502 : 449139 : scalars(1,iffkg)=gxa(1,jj,ia,iproj)
503 : 449139 : scalars(2,iffkg)=gxa(2,jj,ia,iproj)
504 [ + - + - ]: 6048 : else if (choice==2 .and. ndgxdt==1) then
505 : 6048 : jffkg=nffkge+iffkg
506 : : ! Pay attention to the use of reals and imaginary parts here ...
507 : : ! Also, the gxa and dgxdt arrays are switched, in order
508 : : ! to give the correct combination when multiplying ffkg,
509 : : ! see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]]
510 : 6048 : scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
511 : 6048 : scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
512 : 6048 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
513 : 6048 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
514 [ # # ]: 0 : else if (choice==5) then
515 : 0 : jffkg=nffkge+iffkg
516 : : ! The gxa and dgxdt arrays are switched, in order
517 : : ! to give the correct combination when multiplying ffkg,
518 : 0 : scalars(1,jffkg)= gxa(1,jj,ia,iproj)
519 : 0 : scalars(2,jffkg)= gxa(2,jj,ia,iproj)
520 : 0 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
521 : 0 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
522 : : end if
523 : : end do
524 [ + + ]: 147777 : if(choice==3) then
525 [ + + ]: 111888 : do ii=1,ilang3
526 : 103008 : iffkgs=iffkgs+1
527 : 103008 : jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
528 : 103008 : scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
529 : 111888 : scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
530 : : end do
531 : : end if
532 [ + + - + ]: 261552 : if(ilang>=2 .and. choice==5)then
533 [ # # ]: 0 : do ii=1,((ilang-1)*ilang)/2
534 : 0 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
535 : 0 : iffkgk=iffkgk+1
536 : 0 : scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
537 : 0 : scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
538 : : end do
539 : : end if
540 : : end do
541 : : end if
542 : : end do
543 : :
544 : : ! DEBUG
545 : : ! if(choice==5)then
546 : : ! write(std_out,*)' opernl2 : write gxa(:,...) array '
547 : : ! do ii=1,10
548 : : ! write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
549 : : ! end do
550 : : ! write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
551 : : ! do ii=1,10
552 : : ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
553 : : ! end do
554 : : ! write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
555 : : ! do ii=1,4
556 : : ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
557 : : ! end do
558 : : ! end if
559 : :
560 : : ! do iffkg=1,nffkg
561 : : ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
562 : : ! end do
563 : : ! stop
564 : : ! ENDDEBUG
565 : :
566 : : ! ******* Entering the critical part *********************************
567 : :
568 [ + + ]: 7710952 : do ipw=1,nincpw
569 : 7665442 : teffv(1,ipw)=0.0d0
570 : 7710952 : teffv(2,ipw)=0.0d0
571 : : end do
572 [ + + ]: 609753 : do iffkg=1,nffkg
573 [ + + ]: 92985328 : do ipw=1,nincpw
574 : 92375575 : teffv(1,ipw)=teffv(1,ipw)+ffkg(ipw,iffkg)*scalars(1,iffkg)
575 : 92939818 : teffv(2,ipw)=teffv(2,ipw)+ffkg(ipw,iffkg)*scalars(2,iffkg)
576 : : end do
577 : : end do
578 : : ! Multiplication by the complex conjugate of the phase
579 : : ig=ipw1
580 [ + + ]: 7756462 : do ipw=1,nincpw
581 : 7665442 : vect(1,ig)=vect(1,ig)+teffv(1,ipw)*ph3d(1,ig,iaph3d)+teffv(2,ipw)*ph3d(2,ig,iaph3d)
582 : 7665442 : vect(2,ig)=vect(2,ig)+teffv(2,ipw)*ph3d(1,ig,iaph3d)-teffv(1,ipw)*ph3d(2,ig,iaph3d)
583 : 7710952 : ig=ig+1
584 : : end do
585 : :
586 : : ! ******* Leaving the critical part *********************************
587 : :
588 : : ! End loop on atoms
589 : : end do
590 : :
591 : : ! End sign==-1
592 : : else
593 : :
594 : : ! Problem: sign and/or choice do not make sense
595 : : write(message,'(a,2i10,a,a)')&
596 : 0 : & ' Input sign, choice=',sign,choice,ch10,&
597 : 0 : & ' Are not compatible or allowed. '
598 : 0 : ABI_BUG(message)
599 : : end if
600 : :
601 : : ! End loop on blocks of planewaves
602 : : end do
603 : : !!$OMP END DO
604 [ - + ]: 20868 : ABI_FREE(ffkg)
605 [ - + ]: 20868 : ABI_FREE(kpgx)
606 [ - + ]: 20868 : ABI_FREE(parity)
607 [ - + ]: 20868 : ABI_FREE(scalars)
608 [ - + ]: 20868 : ABI_FREE(teffv)
609 : : !!$OMP END PARALLEL
610 : :
611 : :
612 : : !DEBUG
613 : : !if(choice==5)then
614 : : !write(std_out,*)'opernl2 : write vect(2*npw)'
615 : : !do ipw=1,2
616 : : !write(std_out,*)ipw,vect(1:2,ipw)
617 : : !end do
618 : : !write(std_out,*)'opernl2 : write ph3d'
619 : : !do ipw=1,npw
620 : : !write(std_out,*)ipw,ph3d(1:2,ipw,1)
621 : : !end do
622 : : !write(std_out,*)' opernl2 : write gxa array '
623 : : !write(std_out,*)' ang mom ,ia '
624 : : !do iproj=1,mproj
625 : : !do ia=1,1
626 : : !do ii=1,3
627 : : !do ii=1,10
628 : : !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
629 : : !end do
630 : : !end do
631 : : !end do
632 : : !end if
633 : : !if(choice==5)then
634 : : !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
635 : : !do ii=1,10
636 : : !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
637 : : !end do
638 : : !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
639 : : !do ii=1,4
640 : : !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
641 : : !end do
642 : : !stop
643 : : !end if
644 : : !ENDDEBUG
645 : :
646 : 20868 : end subroutine opernl2
647 : : !!***
648 : :
649 : : !!****f* ABINIT/opernl3
650 : : !! NAME
651 : : !! opernl3
652 : : !!
653 : : !! FUNCTION
654 : : !! Operate with the non-local part of the hamiltonian,
655 : : !! either from reciprocal space to projected quantities (sign=1),
656 : : !! or from projected quantities to reciprocal space (sign=-1)
657 : : !!
658 : : !! INPUTS
659 : : !! if(sign==-1 .and. (choice==2 .or choice==4 .or. choice==5))
660 : : !! dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
661 : : !! or with respect to ddk
662 : : !! if(sign==-1 .and. choice==3)
663 : : !! dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
664 : : !! ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
665 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
666 : : !! ia3=gives the number of the first atom in the subset presently treated
667 : : !! idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
668 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
669 : : !! ispinor=1 or 2, gives the spinorial component of ffnl to be used
670 : : !! istwf_k=option parameter that describes the storage of wfs
671 : : !! itypat = type of atom, needed for ffnl
672 : : !! jproj(nlang)=number of projectors for each angular momentum
673 : : !! kg_k(3,npw)=integer coords of planewaves in basis sphere
674 : : !! kpg_k(npw,npkg)= (k+G) components and related data
675 : : !! kpt(3)=real components of k point in terms of recip. translations
676 : : !! lmnmax=max. number of (l,n) components over all type of psps
677 : : !! matblk=dimension of the array ph3d
678 : : !! mincat= maximum increment of atoms
679 : : !! mlang1 = dimensions for dgxdis1
680 : : !! mlang3 = one of the dimensions of the array gxa
681 : : !! mlang4 = dimension for dgxds
682 : : !! mlang5 = dimensions for dgxdis2
683 : : !! mlang6 = dimension for d2gxds2
684 : : !! mproj=maximum dimension for number of projection operators for each
685 : : !! angular momentum for nonlocal pseudopotential
686 : : !! ndgxdt=second dimension of dgxdt
687 : : !! nffnl=third dimension of ffnl
688 : : !! nincat = number of atoms in the subset here treated
689 : : !! nkpg=second size of array kpg_k
690 : : !! nlang = number of angular momenta to be treated = 1 + highest ang. mom.
691 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
692 : : !! npw = number of plane waves in reciprocal space
693 : : !! ntypat = number of type of atoms, dimension needed for ffnl
694 : : !! sign : if 1, go from reciprocal space to projected scalars,
695 : : !! if -1, go from projected scalars to reciprocal space.
696 : : !! if(sign==1),
697 : : !! vect(2*npw)=starting vector in reciprocal space
698 : : !! if(sign==-1)
699 : : !! gxa(2,mlang3,nincat,mproj)=modified projected scalars;
700 : : !! NOTE that metric contractions have already been performed on the
701 : : !! arrays gxa if sign=-1
702 : : !! ph3d(2,npw,matblk)=three-dimensional phase factors
703 : : !!
704 : : !! OUTPUT
705 : : !! if(sign==1)
706 : : !! gxa(2,mlang3,mincat,mproj)= projected scalars
707 : : !! if(sign==1 .and. (choice==2 .or choice==4 .or. choice==5 .or. choice==23))
708 : : !! dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
709 : : !! or with respect to ddk
710 : : !! if(sign==1 .and. (choice==3 .or. choice==23))
711 : : !! dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
712 : : !! if(sign==1 .and. choice==6)
713 : : !! dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars
714 : : !! wrt coord. indexed for internal strain
715 : : !! d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars
716 : : !! wrt strain and coord
717 : : !! d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars
718 : : !! wrt strains
719 : : !! if(sign==-1)
720 : : !! vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.
721 : : !!
722 : : !! NOTES
723 : : !! Operate with the non-local part of the hamiltonian for one type of
724 : : !! atom, and within this given type of atom, for a subset
725 : : !! of at most nincat atoms.
726 : : !!
727 : : !! This routine basically replaces getgla (gxa here is the former gla),
728 : : !! except for the calculation of <G|dVnl/dk|C> or strain gradients.
729 : : !!
730 : : !! SOURCE
731 : :
732 : 252692 : subroutine opernl3(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
733 : 252692 : & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
734 : 252692 : & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
735 : : & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
736 [ - + ]: 252692 : & ntypat,ph3d,sign,vect)
737 : :
738 : : !Arguments ------------------------------------
739 : : !scalars
740 : : integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk
741 : : integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt
742 : : integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat,sign
743 : : !arrays
744 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
745 : : integer,intent(in) :: nloalg(3)
746 : : real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3)
747 : : real(dp),intent(in) :: kpg_k(npw,nkpg),kpt(3),ph3d(2,npw,matblk)
748 : : real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj)
749 : : real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
750 : : real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj),vect(:,:)
751 : : real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj)
752 : : real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj)
753 : : real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj)
754 : :
755 : : !Local variables-------------------------------
756 : : !ntens sets the maximum number of independent tensor components
757 : : !over all allowed angular momenta; need 20 for spdf for tensors
758 : : !up to rank 3; to handle stress tensor, need up to rank 5
759 : : !to handle strain 2DER, need up to rank 7
760 : : !scalars
761 : : integer :: ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2,ilang3
762 : : integer :: ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs,mblkpw
763 : : integer :: mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2,nincpw,nproj
764 : : integer :: ntens
765 : : real(dp) :: ai,ar
766 : : real(dp),parameter :: two_pi2=two_pi**2
767 : : character(len=500) :: message
768 : : !arrays
769 : 252692 : integer,allocatable :: parity(:)
770 : : ! real(dp) :: tsec(2)
771 : 252692 : real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
772 : :
773 : : ! *************************************************************************
774 : :
775 : : !call wrtout(std_out,"in opernl3","COLL")
776 : :
777 : : !mblkpw sets the size of blocks of planewaves
778 : 252692 : mblkpw=NLO_MBLKPW
779 : :
780 : : !Get the actual maximum number of projectors
781 [ + + ]: 1105028 : mmproj=maxval(indlmn(3,:,itypat))
782 : :
783 : : !Initialisation before blocking on the plane waves
784 : :
785 [ + + ]: 252692 : if (sign==1) then
786 : : ! Put projected scalars to zero
787 [ + + + + : 11436560 : gxa(:,:,:,1:mmproj)=0.0d0
+ + + + ]
788 [ + + + + : 2993910 : if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
+ + + + +
+ + + ]
789 [ + + + + : 400046 : if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
+ + + + +
+ ]
790 [ + + ]: 126990 : if (choice==6) then
791 [ + + + + : 13312 : dgxdis(:,:,:,1:mmproj)=0.0d0
+ + + + ]
792 [ + + + + : 28672 : d2gxdis(:,:,:,1:mmproj)=0.0d0
+ + + + ]
793 [ + + + + : 38656 : d2gxds2(:,:,:,1:mmproj)=0.0d0
+ + + + ]
794 : : end if
795 : : end if
796 : :
797 : : !Set up dimension of kpgx and allocate
798 : : !ntens sets the maximum number of independent tensor components
799 : : !over all allowed angular momenta; need 20 for spdf for tensors
800 : : !up to rank 3; to handle stress tensor, need up to rank 5
801 : 252692 : ntens=1
802 [ + + + + ]: 252692 : if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23) ntens=4
803 [ + - + + ]: 252692 : if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
804 [ + - + + : 252692 : if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
+ + ]
805 [ + + + - : 252692 : if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
+ + ]
806 [ + + + - : 252692 : if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
+ + + - ]
807 [ + + - + ]: 252692 : if(choice==6 .and. nlang>=3)ntens=84
808 [ + + - + ]: 252692 : if(choice==6 .and. nlang==4)ntens=120
809 : :
810 : : !Set up second dimension of ffkg array, and allocate
811 : 252692 : nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
812 [ + + ]: 643198 : do ilang=1,nlang
813 : : ! Get the number of projectors for that angular momentum
814 : 390506 : nproj=jproj(ilang)
815 : : ! If there is a non-local part, accumulate the number of vectors needed
816 : : ! The variables ilang below are the number of independent tensors of
817 : : ! various ranks, the variable names being more historical than logical.
818 : : ! ilang2=number of rank ilang-1
819 : : ! ilang3=number of rank ilang+1
820 : : ! ilang4=number of rank ilang
821 : : ! ilang5=number of rank ilang+2
822 : : ! ilang6=number of rank ilang+3
823 [ + - ]: 643198 : if(nproj>0)then
824 : 390506 : ilang2=(ilang*(ilang+1))/2
825 : 390506 : nffkge=nffkge+nproj*ilang2
826 [ - + ]: 390506 : if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
827 [ + + ]: 390506 : if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2
828 [ + + ]: 390506 : if(choice==3 .or. choice==6 .or. choice==23)then
829 : 2960 : ilang3=((ilang+2)*(ilang+3))/2
830 : 2960 : nffkgs=nffkgs+nproj*ilang3
831 : : end if
832 [ + + ]: 390506 : if(choice==6)then
833 : 128 : ilang4=((ilang+1)*(ilang+2))/2
834 : 128 : ilang5=((ilang+3)*(ilang+4))/2
835 : 128 : ilang6=((ilang+4)*(ilang+5))/2
836 : 128 : nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6)
837 : : end if
838 : : end if
839 : : end do
840 : 252692 : nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk
841 : :
842 : : !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a
843 : : !sigfault with nthreads==1
844 : : !
845 : : !Loop on subsets of plane waves (blocking)
846 : : !!$OMP PARALLEL DEFAULT(PRIVATE) &
847 : : !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) &
848 : : !!$OMP SHARED(istwf_k,itypat,jproj,kg_k,kpg_k,kpt,lmnmax,mblkpw,mproj) &
849 : : !!$OMP SHARED(ndgxdt,nffkg,nffkgd,nffkge,nffkgs,nincat,nkpg,nlang) &
850 : : !!$OMP SHARED(nloalg,ph3d,npw,ntens,ntypat,sign,vect)
851 : :
852 [ - + + - : 1010768 : ABI_MALLOC(ffkg,(nffkg,mblkpw))
+ - + - +
- - + -
+ ]
853 [ - + + - : 758076 : ABI_MALLOC(parity,(nffkg))
+ - + - -
+ - + ]
854 [ - + + - : 758076 : ABI_MALLOC(kpgx,(mblkpw,ntens))
+ - + - -
+ - + ]
855 [ - + + - : 758076 : ABI_MALLOC(scalars,(2,nffkg))
+ - + - -
+ - + ]
856 [ - + ]: 252692 : ABI_MALLOC(teffv,(2,mblkpw))
857 : : !!$OMP DO
858 [ + - + - ]: 505384 : do ipw1=1,npw,mblkpw
859 : :
860 : 252692 : ipw2=min(npw,ipw1+mblkpw-1)
861 : 252692 : nincpw=ipw2-ipw1+1
862 : :
863 : : ! Initialize kpgx array related to tensors defined below
864 : : call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
865 : : & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
866 : 252692 : & npw,ntens,ntypat,parity)
867 : :
868 : : ! call timab(74+choice,2,tsec)
869 : :
870 : : ! Now treat the different signs
871 [ + + ]: 505384 : if (sign==1) then
872 : :
873 [ + + ]: 266292 : do ia=1,nincat
874 : :
875 : : ! Compute the shift eventually needed to get the phases in ph3d
876 : 139302 : iaph3d=ia
877 [ + + ]: 139302 : if(nloalg(2)>0)iaph3d=ia+ia3-1
878 : :
879 [ + + ]: 683550 : do iffkg=1,nffkg
880 : 544248 : scalars(1,iffkg)=0.0d0
881 : 683550 : scalars(2,iffkg)=0.0d0
882 : : end do
883 : :
884 : : ! ******* Entering the first time-consuming part of the routine *******
885 : :
886 : : ! Note : first multiply by the phase factor
887 : : ! This allows to be left with only real operations afterwards.
888 : : ig=ipw1
889 : :
890 [ + + ]: 7617946 : do ipw=1,nincpw
891 : 7478644 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
892 : 7478644 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
893 [ + + ]: 43317274 : do iffkg=1,nffkg
894 : 35838630 : scalars(1,iffkg)=scalars(1,iffkg)+ar*ffkg(iffkg,ipw)
895 : 43317274 : scalars(2,iffkg)=scalars(2,iffkg)+ai*ffkg(iffkg,ipw)
896 : : end do
897 : 7617946 : ig=ig+1
898 : : end do
899 : :
900 : : ! ******* Leaving the critical part *********************************
901 : :
902 : : ! DEBUG
903 : : ! write(std_out,*)' opernl2 : scalars'
904 : : ! do iffkg=1,10
905 : : ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
906 : : ! end do
907 : : ! stop
908 : : ! ENDDEBUG
909 : :
910 [ - + ]: 139302 : if(istwf_k>=2)then
911 : : ! Impose parity of resulting scalar (this operation could be
912 : : ! replaced by direct saving of CPU time in the preceeding section)
913 [ # # ]: 0 : do iffkg=1,nffkg
914 : 0 : scalars(parity(iffkg),iffkg)=0.0d0
915 : : end do
916 : : end if
917 : :
918 : 139302 : iffkg=0
919 : 139302 : iffkgs=nffkge+nffkgd
920 : 139302 : iffkgs2=nffkge+nffkgs
921 : 139302 : iffkgk=nffkge*2
922 [ + + ]: 487193 : do ilang=1,nlang
923 : 220901 : nproj=jproj(ilang)
924 [ + - ]: 360203 : if(nproj>0)then
925 : : ! ilang2 is the number of independent tensor components
926 : : ! for symmetric tensor of rank ilang-1
927 : 220901 : ilang2=(ilang*(ilang+1))/2
928 : :
929 : : ! Loop over projectors
930 [ + + ]: 485111 : do iproj=1,nproj
931 : : ! Multiply by the k+G factors (tensors of various rank)
932 [ + + ]: 691618 : do ii=1,ilang2
933 : : ! Get the starting address for the relevant tensor
934 : 427408 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
935 : 427408 : iffkg=iffkg+1
936 : : ! !$OMP CRITICAL (OPERNL3_1)
937 : 427408 : gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
938 : 427408 : gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg)
939 : : ! !$OMP END CRITICAL (OPERNL3_1)
940 : : ! Now, compute gradients, if needed.
941 [ + + + + ]: 427408 : if ((choice==2.or.choice==23) .and. ndgxdt==3) then
942 [ + + ]: 5728 : do mu=1,3
943 : 4296 : jffkg=nffkge+(iffkg-1)*3+mu
944 : : ! Pay attention to the use of reals and imaginary parts here ...
945 : : ! !$OMP CRITICAL (OPERNL3_2)
946 : 4296 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
947 : 5728 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
948 : : ! !$OMP END CRITICAL (OPERNL3_2)
949 : : end do
950 : : end if
951 [ + + + - ]: 427408 : if (choice==2 .and. ndgxdt==1) then
952 : 50832 : jffkg=nffkge+iffkg
953 : : ! Pay attention to the use of reals and imaginary parts here ...
954 : : ! !$OMP CRITICAL (OPERNL3_3)
955 : 50832 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg)
956 : 50832 : dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg)
957 : : ! !$OMP END CRITICAL (OPERNL3_3)
958 : : end if
959 [ + + ]: 691618 : if (choice==4) then
960 [ + + ]: 11264 : do mu=1,3
961 : 8448 : jffkg=nffkge+(iffkg-1)*9+mu
962 : : ! Pay attention to the use of reals and imaginary parts here ...
963 : : ! !$OMP CRITICAL (OPERNL3_4)
964 : 8448 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
965 : 11264 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
966 : : ! !$OMP END CRITICAL (OPERNL3_4)
967 : : end do
968 [ + + ]: 19712 : do mu=4,9
969 : 16896 : jffkg=nffkge+(iffkg-1)*9+mu
970 : : ! Pay attention to the use of reals and imaginary parts here ...
971 : : ! Also, note the multiplication by (2 pi)**2
972 : : ! !$OMP CRITICAL (OPERNL3_5)
973 : 16896 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg)
974 : 19712 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg)
975 : : ! !$OMP END CRITICAL (OPERNL3_5)
976 : : end do
977 : : end if
978 : : ! End loop on ii=1,ilang2
979 : : end do
980 : :
981 [ + + ]: 264210 : if (choice==3 .or. choice==6 .or. choice==23) then
982 : : ! Compute additional tensors related to strain gradients
983 : : ! ilang3 is number of unique tensor components of rank ilang+1
984 : 3432 : ilang3=((ilang+2)*(ilang+3))/2
985 : 3432 : jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
986 : : ! Compute strain gradient tensor components
987 [ + + ]: 30840 : do ii=1,ilang3
988 : : ! Note that iffkgs is also used by ddk and 2nd derivative parts
989 : 27408 : iffkgs=iffkgs+1
990 : 27408 : jj=ii+jjs
991 : : ! !$OMP CRITICAL (OPERNL3_6)
992 : 27408 : dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
993 : 30840 : dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs)
994 : : ! !$OMP END CRITICAL (OPERNL3_6)
995 : : end do
996 : : end if
997 : :
998 [ + + ]: 264210 : if (choice==6) then
999 : : ! Compute additional tensors related to strain 2nd derivatives
1000 : : ! and internal strain derivatives
1001 : : ! ilang6 is number of unique tensor components of rank ilang+3
1002 : 256 : ilang6=((ilang+4)*(ilang+5))/2
1003 : 256 : jjs=((ilang+3)*(ilang+4)*(ilang+5))/6
1004 : : ! Compute strain gradient tensor components
1005 [ + + ]: 4864 : do ii=1,ilang6
1006 : : ! Note that iffkgs is also used by ddk part
1007 : 4608 : iffkgs2=iffkgs2+1
1008 : 4608 : jj=ii+jjs
1009 : : ! !$OMP CRITICAL (OPERNL3_6)
1010 : 4608 : d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2)
1011 : 4864 : d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2)
1012 : : ! !$OMP END CRITICAL (OPERNL3_6)
1013 : : end do
1014 : :
1015 : : ! ilang4 is number of unique tensor components of rank ilang
1016 : 256 : ilang4=((ilang+1)*(ilang+2))/2
1017 : 256 : jjs=((ilang)*(ilang+1)*(ilang+2))/6
1018 : : ! Compute internal strain gradient tensor components
1019 [ + + ]: 1408 : do ii=1,ilang4
1020 : 1152 : iffkgs2=iffkgs2+1
1021 : 1152 : jj=ii+jjs
1022 : : ! !$OMP CRITICAL (OPERNL3_6)
1023 : : ! Pay attention to the use of reals and imaginary parts here ...
1024 : 1152 : dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2)
1025 : 1408 : dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2)
1026 : : ! !$OMP END CRITICAL (OPERNL3_6)
1027 : : end do
1028 : :
1029 : : ! ilang5 is number of unique tensor components of rank ilang+2
1030 : 256 : ilang5=((ilang+3)*(ilang+4))/2
1031 : 256 : jjs=((ilang+2)*(ilang+3)*(ilang+4))/6
1032 : : ! Compute internal strain gradient tensor components
1033 [ + + ]: 3456 : do ii=1,ilang5
1034 : 3200 : iffkgs2=iffkgs2+1
1035 : 3200 : jj=ii+jjs
1036 : : ! !$OMP CRITICAL (OPERNL3_6)
1037 : : ! Pay attention to the use of reals and imaginary parts here ...
1038 : 3200 : d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2)
1039 : 3456 : d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2)
1040 : : ! !$OMP END CRITICAL (OPERNL3_6)
1041 : : end do
1042 : : end if ! choice==6
1043 : :
1044 [ - + ]: 485111 : if (choice==5) then
1045 : : ! Compute additional tensors related to ddk with ffnl(:,2,...)
1046 : 0 : ilangx=(ilang*(ilang+1))/2
1047 : 0 : jjs=((ilang-1)*ilang*(ilang+1))/6
1048 [ # # ]: 0 : do ii=1,ilangx
1049 : : ! Note that iffkgs is also used by strain part
1050 : 0 : iffkgs=iffkgs+1
1051 : 0 : jj=ii+jjs
1052 : : ! !$OMP CRITICAL (OPERNL3_7)
1053 : 0 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
1054 : 0 : dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs)
1055 : : ! !$OMP END CRITICAL (OPERNL3_7)
1056 : : end do
1057 : : ! Compute additional tensors related to ddk with ffnl(:,1,...)
1058 [ # # ]: 0 : if(ilang>=2)then
1059 : 0 : ilangx=((ilang-1)*ilang)/2
1060 : 0 : jjs=((ilang-2)*(ilang-1)*ilang)/6
1061 [ # # ]: 0 : do ii=1,ilangx
1062 : 0 : iffkgk=iffkgk+1
1063 : 0 : jj=ii+jjs
1064 : : ! !$OMP CRITICAL (OPERNL3_8)
1065 : 0 : dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
1066 : 0 : dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk)
1067 : : ! !$OMP END CRITICAL (OPERNL3_8)
1068 : : end do
1069 : : end if
1070 : : end if
1071 : :
1072 : : ! End projector loop
1073 : : end do
1074 : :
1075 : : ! End condition of non-zero projectors
1076 : : end if
1077 : :
1078 : : ! End angular momentum loop
1079 : : end do
1080 : :
1081 : : ! End loop on atoms
1082 : : end do
1083 : :
1084 [ + - ]: 125702 : else if (sign==-1) then
1085 : : ! Application of non-local part from projected scalars
1086 : : ! back to reciprocal space ...
1087 : : ! [this section merely computes terms which add to <G|Vnl|C>;
1088 : : ! nothing here is needed when various gradients are being computed]
1089 : :
1090 : : ! Loop on atoms
1091 [ + + ]: 263484 : do ia=1,nincat
1092 : :
1093 : : ! Compute the shift eventually needed to get the phases in ph3d
1094 : 137782 : iaph3d=ia
1095 [ + + ]: 137782 : if(nloalg(2)>0)iaph3d=ia+ia3-1
1096 : :
1097 : : ! Transfer gxa (and eventually dgxdt) in scalars with different indexing
1098 : 137782 : iffkg=0
1099 : 137782 : iffkgk=nffkge*2
1100 : 137782 : iffkgs=nffkge
1101 [ + + ]: 356171 : do ilang=1,nlang
1102 : 218389 : nproj=jproj(ilang)
1103 [ + - ]: 356171 : if (nproj>0) then
1104 : 218389 : ilang2=(ilang*(ilang+1))/2
1105 : 218389 : ilang3=((ilang+2)*(ilang+3))/2
1106 [ + + ]: 479823 : do iproj=1,nproj
1107 [ + + ]: 684082 : do ii=1,ilang2
1108 : 422648 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
1109 : 422648 : iffkg=iffkg+1
1110 [ + + ]: 684082 : if(choice==1 .or. choice==3)then
1111 : 371816 : scalars(1,iffkg)=gxa(1,jj,ia,iproj)
1112 : 371816 : scalars(2,iffkg)=gxa(2,jj,ia,iproj)
1113 [ + - + - ]: 50832 : else if (choice==2 .and. ndgxdt==1) then
1114 : 50832 : jffkg=nffkge+iffkg
1115 : : ! Pay attention to the use of reals and imaginary parts here ...
1116 : : ! Also, the gxa and dgxdt arrays are switched, in order
1117 : : ! to give the correct combination when multiplying ffkg,
1118 : : ! see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]]
1119 : 50832 : scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
1120 : 50832 : scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
1121 : 50832 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
1122 : 50832 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
1123 [ # # ]: 0 : else if (choice==5) then
1124 : 0 : jffkg=nffkge+iffkg
1125 : : ! The gxa and dgxdt arrays are switched, in order
1126 : : ! to give the correct combination when multiplying ffkg,
1127 : 0 : scalars(1,jffkg)= gxa(1,jj,ia,iproj)
1128 : 0 : scalars(2,jffkg)= gxa(2,jj,ia,iproj)
1129 : 0 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
1130 : 0 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
1131 : : end if
1132 : : end do
1133 [ + + ]: 261434 : if(choice==3) then
1134 [ + + ]: 22032 : do ii=1,ilang3
1135 : 19584 : iffkgs=iffkgs+1
1136 : 19584 : jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
1137 : 19584 : scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
1138 : 22032 : scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
1139 : : end do
1140 : : end if
1141 [ + + - + ]: 479823 : if(ilang>=2 .and. choice==5)then
1142 [ # # ]: 0 : do ii=1,((ilang-1)*ilang)/2
1143 : 0 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
1144 : 0 : iffkgk=iffkgk+1
1145 : 0 : scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
1146 : 0 : scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
1147 : : end do
1148 : : end if
1149 : : end do
1150 : : end if
1151 : : end do
1152 : :
1153 : : ! DEBUG
1154 : : ! if(choice==5)then
1155 : : ! write(std_out,*)' opernl2 : write gxa(:,...) array '
1156 : : ! do ii=1,10
1157 : : ! write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
1158 : : ! end do
1159 : : ! write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
1160 : : ! do ii=1,10
1161 : : ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
1162 : : ! end do
1163 : : ! write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
1164 : : ! do ii=1,4
1165 : : ! write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
1166 : : ! end do
1167 : : ! end if
1168 : :
1169 : : ! do iffkg=1,nffkg
1170 : : ! write(std_out,*)'iffkg,scalars',iffkg,scalars(1:2,iffkg)
1171 : : ! end do
1172 : : ! stop
1173 : : ! ENDDEBUG
1174 : :
1175 : : ig=ipw1
1176 : :
1177 : : ! ******* Entering the critical part *********************************
1178 : :
1179 [ + + ]: 7620732 : do ipw=1,nincpw
1180 : : ar=0.0d0
1181 : : ai=0.0d0
1182 [ + + ]: 38079142 : do iffkg=1,nffkg
1183 : 30721894 : ar=ar+ffkg(iffkg,ipw)*scalars(1,iffkg)
1184 : 38079142 : ai=ai+ffkg(iffkg,ipw)*scalars(2,iffkg)
1185 : : end do
1186 : 7357248 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
1187 : 7357248 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
1188 : 7495030 : ig=ig+1
1189 : : end do
1190 : :
1191 : : ! ******* Leaving the critical part *********************************
1192 : :
1193 : : ! End loop on atoms
1194 : : end do
1195 : :
1196 : : ! End sign==-1
1197 : : else
1198 : :
1199 : : ! Problem: sign and/or choice do not make sense
1200 : : write(message, '(a,2i10,a,a)' )&
1201 : 0 : & ' Input sign, choice=',sign,choice,ch10,&
1202 : 0 : & ' Are not compatible or allowed. '
1203 : 0 : ABI_BUG(message)
1204 : : end if
1205 : :
1206 : : ! End loop on blocks of planewaves
1207 : : end do
1208 : : !!$OMP END DO
1209 [ - + ]: 252692 : ABI_FREE(ffkg)
1210 [ - + ]: 252692 : ABI_FREE(kpgx)
1211 [ - + ]: 252692 : ABI_FREE(parity)
1212 [ - + ]: 252692 : ABI_FREE(scalars)
1213 [ - + ]: 252692 : ABI_FREE(teffv)
1214 : : !!$OMP END PARALLEL
1215 : :
1216 : :
1217 : : !DEBUG
1218 : : !if(choice==5)then
1219 : : !write(std_out,*)'opernl2 : write vect(2*npw)'
1220 : : !do ipw=1,2
1221 : : !write(std_out,*)ipw,vect(1:2,ipw)
1222 : : !end do
1223 : : !write(std_out,*)'opernl2 : write ph3d'
1224 : : !do ipw=1,npw
1225 : : !write(std_out,*)ipw,ph3d(1:2,ipw,1)
1226 : : !end do
1227 : : !write(std_out,*)' opernl2 : write gxa array '
1228 : : !write(std_out,*)' ang mom ,ia '
1229 : : !do iproj=1,mproj
1230 : : !do ia=1,1
1231 : : !do ii=1,3
1232 : : !do ii=1,10
1233 : : !write(std_out,'(i3,2es16.6)' )ii,gxa(:,ii,1,1)
1234 : : !end do
1235 : : !end do
1236 : : !end do
1237 : : !end if
1238 : : !if(choice==5)then
1239 : : !write(std_out,*)' opernl2 : write dgxdt(:,1,...) array '
1240 : : !do ii=1,10
1241 : : !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,1,ii,1,1)
1242 : : !end do
1243 : : !write(std_out,*)' opernl2 : write dgxdt(:,2,...) array '
1244 : : !do ii=1,4
1245 : : !write(std_out,'(i3,2es16.6)' )ii,dgxdt(:,2,ii,1,1)
1246 : : !end do
1247 : : !stop
1248 : : !end if
1249 : : !ENDDEBUG
1250 : :
1251 : 252692 : end subroutine opernl3
1252 : : !!***
1253 : :
1254 : :
1255 : : !!****f* ABINIT/opernl4a
1256 : : !! NAME
1257 : : !! opernl4a
1258 : : !!
1259 : : !! FUNCTION
1260 : : !! Operate with the non-local part of the hamiltonian,
1261 : : !! from reciprocal space to projected quantities
1262 : : !! (oprnl4b is from projected quantities to reciprocal space)
1263 : : !!
1264 : : !! INPUTS
1265 : : !! ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
1266 : : !! ia3=gives the number of the first atom in the subset presently treated
1267 : : !! idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
1268 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
1269 : : !! ispinor=1 or 2, gives the spinorial component of ffnl to be used
1270 : : !! istwf_k=option parameter that describes the storage of wfs
1271 : : !! itypat = type of atom, needed for ffnl
1272 : : !! jproj(nlang)=number of projectors for each angular momentum
1273 : : !! kg_k(3,npw)=integer coords of planewaves in basis sphere
1274 : : !! kpg_k(npw,npkg)= (k+G) components and related data
1275 : : !! kpt(3)=real components of k point in terms of recip. translations
1276 : : !! lmnmax=max. number of (l,n) components over all type of psps
1277 : : !! matblk=dimension of the array ph3d
1278 : : !! mincat= maximum increment of atoms
1279 : : !! mlang1 = dimensions for dgxdis1
1280 : : !! mlang3 = one of the dimensions of the array gxa
1281 : : !! mlang4 = dimension for dgxds
1282 : : !! mlang5 = dimensions for dgxdis2
1283 : : !! mlang6 = dimension for d2gxds2
1284 : : !! mproj=maximum dimension for number of projection operators for each
1285 : : !! angular momentum for nonlocal pseudopotential
1286 : : !! ndgxdt=second dimension of dgxdt
1287 : : !! nffnl=third dimension of ffnl
1288 : : !! nincat = number of atoms in the subset here treated
1289 : : !! nkpg=second size of array kpg_k
1290 : : !! nlang = number of angular momenta to be treated = 1 + highest ang. mom.
1291 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
1292 : : !! npw = number of plane waves in reciprocal space
1293 : : !! ntypat = number of type of atoms, dimension needed for ffnl
1294 : : !! vect(2*npw)=starting vector in reciprocal space
1295 : : !! ph3d(2,npw,matblk)=three-dimensional phase factors
1296 : : !!
1297 : : !! OUTPUT
1298 : : !! gxa(2,mlang3,mincat,mproj)= projected scalars
1299 : : !! if(choice==2 .or choice==4 .or. choice==5 .or. choice==23)
1300 : : !! dgxdt(2,ndgxdt,mlang3,mincat,mproj)= gradients of projected scalars wrt coords
1301 : : !! or with respect to ddk
1302 : : !! if(choice==3 .or. choice==23)
1303 : : !! dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
1304 : : !! if(choice==6)
1305 : : !! dgxdis((2,mlang1,mincat,mproj) = derivatives of projected scalars
1306 : : !! wrt coord. indexed for internal strain
1307 : : !! d2gxdis((2,mlang5,mincat,mproj) = 2nd derivatives of projected scalars
1308 : : !! wrt strain and coord
1309 : : !! d2gxds2((2,mlang6,mincat,mproj) = 2nd derivatives of projected scalars
1310 : : !! wrt strains
1311 : : !!
1312 : : !! NOTES
1313 : : !! Operate with the non-local part of the hamiltonian for one type of
1314 : : !! atom, and within this given type of atom, for a subset
1315 : : !! of at most nincat atoms.
1316 : : !! This routine basically replaces getgla (gxa here is the former gla),
1317 : : !! except for the calculation of <G|dVnl/dk|C> or strain gradients.
1318 : : !!
1319 : : !! Present version decomposed according to iffkg
1320 : : !!
1321 : : !! SOURCE
1322 : :
1323 : 48511395 : subroutine opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
1324 : 48511395 : & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
1325 : 48511395 : & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1326 : : & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
1327 [ - + ]: 48511395 : & ntypat,ph3d,vect)
1328 : :
1329 : : !Arguments ------------------------------------
1330 : : !scalars
1331 : : integer,intent(in) :: choice,ia3,idir,ispinor,istwf_k,itypat,lmnmax,matblk
1332 : : integer,intent(in) :: mincat,mlang1,mlang3,mlang4,mlang5,mlang6,mproj,ndgxdt
1333 : : integer,intent(in) :: nffnl,nincat,nkpg,nlang,npw,ntypat
1334 : : !arrays
1335 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
1336 : : integer,intent(in) :: nloalg(3)
1337 : : real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
1338 : : real(dp),intent(in) :: kpt(3),ph3d(2,npw,matblk),vect(:,:)
1339 : : real(dp),intent(out) :: d2gxdis(2,mlang5,mincat,mproj)
1340 : : real(dp),intent(out) :: d2gxds2(2,mlang6,mincat,mproj)
1341 : : real(dp),intent(out) :: dgxdis(2,mlang1,mincat,mproj)
1342 : : real(dp),intent(inout) :: dgxds(2,mlang4,mincat,mproj) !vz_i
1343 : : real(dp),intent(inout) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj) !vz_i
1344 : : real(dp),intent(inout) :: gxa(2,mlang3,mincat,mproj) !vz_i
1345 : :
1346 : : !Local variables-------------------------------
1347 : : !scalars
1348 : : integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,iffkgs2,ig,ii,ilang,ilang2
1349 : : integer :: ilang3,ilang4,ilang5,ilang6,ilangx,iproj,ipw,ipw1,ipw2,jffkg,jj,jjs
1350 : : integer :: jump,mblkpw,mmproj,mu,nffkg,nffkgd,nffkge,nffkgk,nffkgs,nffkgs2
1351 : : integer :: nincpw,nproj,ntens,start
1352 : : real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8
1353 : : real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8
1354 : : real(dp),parameter :: two_pi2=two_pi*two_pi
1355 : : !arrays
1356 : 48511395 : integer,allocatable :: parity(:)
1357 : 48511395 : real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
1358 : :
1359 : : ! *************************************************************************
1360 : :
1361 : : !call wrtout(std_out,"in opernl4a","COLL")
1362 : :
1363 : : !mblkpw sets the size of blocks of planewaves
1364 : 48511395 : mblkpw=NLO_MBLKPW
1365 : :
1366 : : !jump governs, in fine, the use of registers in the most cpu
1367 : : !time consuming part of the routine. Until now, jump=8 is the maximal value.
1368 : : !The optimal value will be machine-dependent !
1369 : 48511395 : jump=4
1370 : :
1371 : : !Get the actual maximum number of projectors
1372 [ + + ]: 343907992 : mmproj=maxval(indlmn(3,:,itypat))
1373 : :
1374 : : !Initialisation before blocking on the plane waves
1375 : :
1376 : : !Put projected scalars to zero
1377 [ + + + + : 10463140655 : gxa(:,:,:,1:mmproj)=0.0d0
+ + + + ]
1378 [ + + + + : 1895517987 : if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
+ + + + +
+ + + ]
1379 [ + + + + : 496149777 : if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
+ + + + +
+ ]
1380 [ + + ]: 48511395 : if (choice==6) then
1381 [ + + + + : 5914160 : dgxdis(:,:,:,1:mmproj)=0.0d0
+ + + + ]
1382 [ + + + + : 12709040 : d2gxdis(:,:,:,1:mmproj)=0.0d0
+ + + + ]
1383 [ + + + + : 17125712 : d2gxds2(:,:,:,1:mmproj)=0.0d0
+ + + + ]
1384 : : end if
1385 : :
1386 : : !Set up dimension of kpgx and allocate
1387 : : !ntens sets the maximum number of independent tensor components
1388 : : !over all allowed angular momenta; need 20 for spdf for tensors
1389 : : !up to rank 3; to handle stress tensor, need up to rank 5
1390 : 48511395 : ntens=1
1391 [ + + ]: 49735480 : if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23)ntens=4
1392 [ + + + + ]: 48511395 : if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
1393 [ + + + + : 48511395 : if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
+ + ]
1394 [ + + + + : 48511395 : if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
+ + ]
1395 [ + + + + : 48511395 : if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
+ + + + ]
1396 [ + + + + ]: 48511395 : if(choice==6 .and. nlang>=3)ntens=84
1397 [ + + - + ]: 48511395 : if(choice==6 .and. nlang==4)ntens=120
1398 : :
1399 : : !Set up second dimension of ffkg array, and allocate
1400 : 48511395 : nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
1401 [ + + ]: 167187413 : do ilang=1,nlang
1402 : : ! Get the number of projectors for that angular momentum
1403 : 118676018 : nproj=jproj(ilang)
1404 : : ! If there is a non-local part, accumulate the number of vectors needed
1405 : : ! The variables ilang below are the number of independent tensors of
1406 : : ! various ranks, the variable names being more historical than logical.
1407 : : ! ilang2=number of rank ilang-1
1408 : : ! ilang3=number of rank ilang+1
1409 : : ! ilang4=number of rank ilang
1410 : : ! ilang5=number of rank ilang+2
1411 : : ! ilang6=number of rank ilang+3
1412 [ + + ]: 167187413 : if(nproj>0)then
1413 : 112159501 : ilang2=(ilang*(ilang+1))/2
1414 : 112159501 : nffkge=nffkge+nproj*ilang2
1415 [ + + ]: 112159501 : if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
1416 : 13952727 : if(choice==2 .or. choice==4 .or. choice==23)nffkgd=nffkgd+ndgxdt*nproj*ilang2
1417 : : if(choice==3 .or. choice==6 .or. choice==23)then
1418 : 2887295 : ilang3=((ilang+2)*(ilang+3))/2
1419 : 2887295 : nffkgs=nffkgs+nproj*ilang3
1420 : : end if
1421 [ + + ]: 112159501 : if(choice==6)then
1422 : 80200 : ilang4=((ilang+1)*(ilang+2))/2
1423 : 80200 : ilang5=((ilang+3)*(ilang+4))/2
1424 : 80200 : ilang6=((ilang+4)*(ilang+5))/2
1425 : 80200 : nffkgs2=nffkgs2+nproj*(ilang4+ilang5+ilang6)
1426 : : end if
1427 : : end if
1428 : : end do
1429 : 48511395 : nffkg=nffkge+nffkgd+nffkgs+nffkgs2+nffkgk
1430 : :
1431 : : !DEBUG
1432 : : !write(std_out,*)' jproj(1:nlang)',jproj(1:nlang)
1433 : : !write(std_out,*)' nffkg,nffkge,nffkgd,nffkgs,nffkgk',nffkg,nffkge,nffkgd,nffkgs,nffkgk
1434 : : !ENDDEBUG
1435 : :
1436 : : !Loop on subsets of plane waves (blocking)
1437 : :
1438 : : !Disabled by MG on Dec 6 2011, omp sections have to be tested, this coding causes a
1439 : : !sigfault with nthreads==1
1440 : : !Feb 16 2012: The code does not crash anymore but it's not efficient.
1441 : : !
1442 : : !!$OMP PARALLEL DEFAULT(PRIVATE) &
1443 : : !!$OMP SHARED (choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt) &
1444 : : !!$OMP SHARED (ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat) &
1445 : : !!$OMP SHARED (jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4) &
1446 : : !!$OMP SHARED (mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw) &
1447 : : !!$OMP SHARED (ntypat,ph3d,vect) &
1448 : : !!$OMP SHARED (mblkpw,jump,nffkgd,nffkg,nffkge,nffkgs,ntens)
1449 : :
1450 [ - + + - : 194045580 : ABI_MALLOC(ffkg,(nffkg,mblkpw))
+ - + - +
- - + -
+ ]
1451 [ - + + - : 145534185 : ABI_MALLOC(parity,(nffkg))
+ - + - -
+ - + ]
1452 [ - + + - : 145534185 : ABI_MALLOC(kpgx,(mblkpw,ntens))
+ - + - -
+ - + ]
1453 [ - + + - : 145534185 : ABI_MALLOC(scalars,(2,nffkg))
+ - + - -
+ - + ]
1454 [ - + ]: 48511395 : ABI_MALLOC(teffv,(2,mblkpw))
1455 : :
1456 : : !!$OMP DO
1457 [ + - + + ]: 141649771 : do ipw1=1,npw,mblkpw
1458 : :
1459 : 93138376 : ipw2=min(npw,ipw1+mblkpw-1)
1460 : 93138376 : nincpw=ipw2-ipw1+1
1461 : :
1462 : : ! Initialize kpgx array related to tensors defined below
1463 : : call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
1464 : : & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
1465 : 93138376 : & npw,ntens,ntypat,parity)
1466 : :
1467 [ + + ]: 288957783 : do ia=1,nincat
1468 : :
1469 : : ! Compute the shift eventually needed to get the phases in ph3d
1470 : 147308012 : iaph3d=ia
1471 [ + + ]: 147308012 : if(nloalg(2)>0)iaph3d=ia+ia3-1
1472 : :
1473 [ + + ]: 2531143581 : do iffkg=1,nffkg
1474 : 2531143581 : scalars(1,iffkg)=0.0d0 ; scalars(2,iffkg)=0.0d0
1475 : : end do
1476 : :
1477 : : ! DEBUG
1478 : : ! write(std_out,*)'opernl4, before first time-consuming'
1479 : : ! write(std_out,*)'opernl4 : nffkg,nincpw=',nffkg,nincpw
1480 : : ! write(std_out,*)'ig,ipw,ffkg(1:4),vect(1:2)'
1481 : : ! ig=ipw1
1482 : : ! do ipw=1,nincpw
1483 : : ! write(std_out,'(2i4,13es11.3)' )ig,ipw,ffkg(1:min(9,nffkg),ipw),vect(1:2,ipw),ph3d(1:2,ipw,iaph3d)
1484 : : ! ig=ig+1
1485 : : ! end do
1486 : : ! stop
1487 : : ! ENDDEBUG
1488 : :
1489 : : ! ******* Entering the first time-consuming part of the routine *******
1490 : :
1491 : :
1492 : : ! First, treat small nffkg; send treat the initial phase of big
1493 : : ! nffkg; finally treat the loop needed for big nffkg
1494 : :
1495 : : ! In the loops, first multiply by the phase factor.
1496 : : ! This allows to be left with only real operations afterwards.
1497 : :
1498 : : ! For the time being, the maximal jump allowed is 8.
1499 : :
1500 : : ! 1) Here, treat small nffkg
1501 [ + + ]: 147308012 : if(nffkg<=jump)then
1502 : :
1503 : : select case(nffkg)
1504 : :
1505 : : case(1)
1506 : :
1507 : : scr1=0.0d0 ; sci1=0.0d0
1508 : : ig=ipw1
1509 [ + + ]: 1792229306 : do ipw=1,nincpw
1510 : 1779167587 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1511 : 1779167587 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1512 : 1779167587 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1513 : 1792229306 : ig=ig+1
1514 : : end do
1515 : 13061719 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1516 : :
1517 : : case(2)
1518 : :
1519 : : ig=ipw1
1520 : : scr1=0.0d0 ; sci1=0.0d0
1521 : : scr2=0.0d0 ; sci2=0.0d0
1522 [ + + ]: 320185962 : do ipw=1,nincpw
1523 : 317823230 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1524 : 317823230 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1525 : 317823230 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1526 : 317823230 : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1527 : 320185962 : ig=ig+1
1528 : : end do
1529 : 2362732 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1530 : 2362732 : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1531 : :
1532 : : case(3)
1533 : :
1534 : : ig=ipw1
1535 : : scr1=0.0d0 ; sci1=0.0d0
1536 : : scr2=0.0d0 ; sci2=0.0d0
1537 : : scr3=0.0d0 ; sci3=0.0d0
1538 [ + + ]: 1359388 : do ipw=1,nincpw
1539 : 1352532 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1540 : 1352532 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1541 : 1352532 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1542 : 1352532 : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1543 : 1352532 : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1544 : 1359388 : ig=ig+1
1545 : : end do
1546 : 6856 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1547 : 6856 : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1548 : 6856 : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1549 : :
1550 : : case(4)
1551 : :
1552 : : ig=ipw1
1553 : : scr1=0.0d0 ; sci1=0.0d0
1554 : : scr2=0.0d0 ; sci2=0.0d0
1555 : : scr3=0.0d0 ; sci3=0.0d0
1556 : : scr4=0.0d0 ; sci4=0.0d0
1557 [ + + ]: 1580670979 : do ipw=1,nincpw
1558 : 1567963204 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1559 : 1567963204 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1560 : 1567963204 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1561 : 1567963204 : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1562 : 1567963204 : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1563 : 1567963204 : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1564 : 1580670979 : ig=ig+1
1565 : : end do
1566 : 12707775 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1567 : 12707775 : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1568 : 12707775 : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1569 : 12707775 : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1570 : :
1571 : : case(5)
1572 : :
1573 : : ig=ipw1
1574 : : scr1=0.0d0 ; sci1=0.0d0
1575 : : scr2=0.0d0 ; sci2=0.0d0
1576 : : scr3=0.0d0 ; sci3=0.0d0
1577 : : scr4=0.0d0 ; sci4=0.0d0
1578 : : scr5=0.0d0 ; sci5=0.0d0
1579 : : do ipw=1,nincpw
1580 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1581 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1582 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1583 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1584 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1585 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1586 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1587 : : ig=ig+1
1588 : : end do
1589 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1590 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1591 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1592 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1593 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1594 : :
1595 : : case(6)
1596 : :
1597 : : ig=ipw1
1598 : : scr1=0.0d0 ; sci1=0.0d0
1599 : : scr2=0.0d0 ; sci2=0.0d0
1600 : : scr3=0.0d0 ; sci3=0.0d0
1601 : : scr4=0.0d0 ; sci4=0.0d0
1602 : : scr5=0.0d0 ; sci5=0.0d0
1603 : : scr6=0.0d0 ; sci6=0.0d0
1604 : : do ipw=1,nincpw
1605 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1606 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1607 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1608 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1609 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1610 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1611 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1612 : : scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1613 : : ig=ig+1
1614 : : end do
1615 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1616 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1617 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1618 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1619 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1620 : : scalars(1,6)=scr6 ; scalars(2,6)=sci6
1621 : :
1622 : : case(7)
1623 : :
1624 : : ig=ipw1
1625 : : scr1=0.0d0 ; sci1=0.0d0
1626 : : scr2=0.0d0 ; sci2=0.0d0
1627 : : scr3=0.0d0 ; sci3=0.0d0
1628 : : scr4=0.0d0 ; sci4=0.0d0
1629 : : scr5=0.0d0 ; sci5=0.0d0
1630 : : scr6=0.0d0 ; sci6=0.0d0
1631 : : scr7=0.0d0 ; sci7=0.0d0
1632 : : do ipw=1,nincpw
1633 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1634 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1635 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1636 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1637 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1638 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1639 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1640 : : scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1641 : : scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1642 : : ig=ig+1
1643 : : end do
1644 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1645 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1646 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1647 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1648 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1649 : : scalars(1,6)=scr6 ; scalars(2,6)=sci6
1650 : : scalars(1,7)=scr7 ; scalars(2,7)=sci7
1651 : :
1652 : : case(8)
1653 : :
1654 : : ig=ipw1
1655 : : scr1=0.0d0 ; sci1=0.0d0
1656 : : scr2=0.0d0 ; sci2=0.0d0
1657 : : scr3=0.0d0 ; sci3=0.0d0
1658 : : scr4=0.0d0 ; sci4=0.0d0
1659 : : scr5=0.0d0 ; sci5=0.0d0
1660 : : scr6=0.0d0 ; sci6=0.0d0
1661 : : scr7=0.0d0 ; sci7=0.0d0
1662 : : scr8=0.0d0 ; sci8=0.0d0
1663 : : do ipw=1,nincpw
1664 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1665 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1666 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1667 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1668 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1669 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1670 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1671 : : scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1672 : : scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1673 : : scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw)
1674 : : ig=ig+1
1675 : : end do
1676 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1677 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1678 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1679 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1680 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1681 : : scalars(1,6)=scr6 ; scalars(2,6)=sci6
1682 : : scalars(1,7)=scr7 ; scalars(2,7)=sci7
1683 [ + + + + : 28139082 : scalars(1,8)=scr8 ; scalars(2,8)=sci8
- ]
1684 : :
1685 : : end select
1686 : :
1687 : : else
1688 : : ! Now treat big nffkg
1689 : :
1690 : : ! 2) Here, initialize big nffkg. The only difference with the
1691 : : ! preceeding case is that the intermediate results are stored.
1692 : :
1693 : : select case(jump)
1694 : :
1695 : : case(1)
1696 : :
1697 : : scr1=0.0d0 ; sci1=0.0d0
1698 : : ig=ipw1
1699 : : do ipw=1,nincpw
1700 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1701 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1702 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1703 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1704 : : ig=ig+1
1705 : : end do
1706 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1707 : :
1708 : : case(2)
1709 : :
1710 : : ig=ipw1
1711 : : scr1=0.0d0 ; sci1=0.0d0
1712 : : scr2=0.0d0 ; sci2=0.0d0
1713 : : do ipw=1,nincpw
1714 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1715 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1716 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1717 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1718 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1719 : : ig=ig+1
1720 : : end do
1721 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1722 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1723 : :
1724 : : case(3)
1725 : :
1726 : : ig=ipw1
1727 : : scr1=0.0d0 ; sci1=0.0d0
1728 : : scr2=0.0d0 ; sci2=0.0d0
1729 : : scr3=0.0d0 ; sci3=0.0d0
1730 : : do ipw=1,nincpw
1731 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1732 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1733 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1734 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1735 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1736 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1737 : : ig=ig+1
1738 : : end do
1739 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1740 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1741 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1742 : :
1743 : : case(4)
1744 : :
1745 : : ig=ipw1
1746 : : scr1=0.0d0 ; sci1=0.0d0
1747 : : scr2=0.0d0 ; sci2=0.0d0
1748 : : scr3=0.0d0 ; sci3=0.0d0
1749 : : scr4=0.0d0 ; sci4=0.0d0
1750 [ + + ]: 18648797052 : do ipw=1,nincpw
1751 : 18529628122 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1752 : 18529628122 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1753 : 18529628122 : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1754 : 18529628122 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1755 : 18529628122 : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1756 : 18529628122 : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1757 : 18529628122 : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1758 : 18648797052 : ig=ig+1
1759 : : end do
1760 : 119168930 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1761 : 119168930 : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1762 : 119168930 : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1763 : 119168930 : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1764 : :
1765 : : case(5)
1766 : :
1767 : : ig=ipw1
1768 : : scr1=0.0d0 ; sci1=0.0d0
1769 : : scr2=0.0d0 ; sci2=0.0d0
1770 : : scr3=0.0d0 ; sci3=0.0d0
1771 : : scr4=0.0d0 ; sci4=0.0d0
1772 : : scr5=0.0d0 ; sci5=0.0d0
1773 : : do ipw=1,nincpw
1774 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1775 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1776 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1777 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1778 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1779 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1780 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1781 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1782 : : ig=ig+1
1783 : : end do
1784 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1785 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1786 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1787 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1788 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1789 : :
1790 : : case(6)
1791 : :
1792 : : ig=ipw1
1793 : : scr1=0.0d0 ; sci1=0.0d0
1794 : : scr2=0.0d0 ; sci2=0.0d0
1795 : : scr3=0.0d0 ; sci3=0.0d0
1796 : : scr4=0.0d0 ; sci4=0.0d0
1797 : : scr5=0.0d0 ; sci5=0.0d0
1798 : : scr6=0.0d0 ; sci6=0.0d0
1799 : : do ipw=1,nincpw
1800 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1801 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1802 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1803 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1804 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1805 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1806 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1807 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1808 : : scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1809 : : ig=ig+1
1810 : : end do
1811 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1812 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1813 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1814 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1815 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1816 : : scalars(1,6)=scr6 ; scalars(2,6)=sci6
1817 : :
1818 : : case(7)
1819 : :
1820 : : ig=ipw1
1821 : : scr1=0.0d0 ; sci1=0.0d0
1822 : : scr2=0.0d0 ; sci2=0.0d0
1823 : : scr3=0.0d0 ; sci3=0.0d0
1824 : : scr4=0.0d0 ; sci4=0.0d0
1825 : : scr5=0.0d0 ; sci5=0.0d0
1826 : : scr6=0.0d0 ; sci6=0.0d0
1827 : : scr7=0.0d0 ; sci7=0.0d0
1828 : : do ipw=1,nincpw
1829 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1830 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1831 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1832 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1833 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1834 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1835 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1836 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1837 : : scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1838 : : scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1839 : : ig=ig+1
1840 : : end do
1841 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1842 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1843 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1844 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1845 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1846 : : scalars(1,6)=scr6 ; scalars(2,6)=sci6
1847 : : scalars(1,7)=scr7 ; scalars(2,7)=sci7
1848 : :
1849 : : case(8)
1850 : :
1851 : : ig=ipw1
1852 : : scr1=0.0d0 ; sci1=0.0d0
1853 : : scr2=0.0d0 ; sci2=0.0d0
1854 : : scr3=0.0d0 ; sci3=0.0d0
1855 : : scr4=0.0d0 ; sci4=0.0d0
1856 : : scr5=0.0d0 ; sci5=0.0d0
1857 : : scr6=0.0d0 ; sci6=0.0d0
1858 : : scr7=0.0d0 ; sci7=0.0d0
1859 : : scr8=0.0d0 ; sci8=0.0d0
1860 : : do ipw=1,nincpw
1861 : : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1862 : : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1863 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1864 : : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1865 : : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1866 : : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1867 : : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1868 : : scr5=scr5+ar*ffkg(5,ipw) ; sci5=sci5+ai*ffkg(5,ipw)
1869 : : scr6=scr6+ar*ffkg(6,ipw) ; sci6=sci6+ai*ffkg(6,ipw)
1870 : : scr7=scr7+ar*ffkg(7,ipw) ; sci7=sci7+ai*ffkg(7,ipw)
1871 : : scr8=scr8+ar*ffkg(8,ipw) ; sci8=sci8+ai*ffkg(8,ipw)
1872 : : ig=ig+1
1873 : : end do
1874 : : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1875 : : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1876 : : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1877 : : scalars(1,4)=scr4 ; scalars(2,4)=sci4
1878 : : scalars(1,5)=scr5 ; scalars(2,5)=sci5
1879 : : scalars(1,6)=scr6 ; scalars(2,6)=sci6
1880 : : scalars(1,7)=scr7 ; scalars(2,7)=sci7
1881 : : scalars(1,8)=scr8 ; scalars(2,8)=sci8
1882 : :
1883 : : end select
1884 : :
1885 : : ! 3) Here, do-loop for big nffkg.
1886 : :
1887 [ + - + + ]: 602310362 : do start=1+jump,nffkg,jump
1888 : 483141432 : chunk=min(jump,nffkg-start+1)
1889 : :
1890 : 119168930 : select case(chunk)
1891 : :
1892 : : case(1)
1893 : :
1894 : : scr1=0.0d0 ; sci1=0.0d0
1895 [ + + ]: 2145854822 : do ipw=1,nincpw
1896 : 2130816252 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1897 : 2145854822 : scr1=scr1+ar*ffkg(start,ipw) ; sci1=sci1+ai*ffkg(start,ipw)
1898 : : end do
1899 : 15038570 : scalars(1,start)=scr1 ; scalars(2,start)=sci1
1900 : :
1901 : : case(2)
1902 : :
1903 : : scr1=0.0d0 ; sci1=0.0d0
1904 : : scr2=0.0d0 ; sci2=0.0d0
1905 [ + + ]: 3642953535 : do ipw=1,nincpw
1906 : 3620374239 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1907 : 3620374239 : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1908 : 3642953535 : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1909 : : end do
1910 : 22579296 : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1911 : 22579296 : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1912 : :
1913 : : case(3)
1914 : :
1915 : : scr1=0.0d0 ; sci1=0.0d0
1916 : : scr2=0.0d0 ; sci2=0.0d0
1917 : : scr3=0.0d0 ; sci3=0.0d0
1918 [ + + ]: 590381833 : do ipw=1,nincpw
1919 : 586611405 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1920 : 586611405 : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1921 : 586611405 : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1922 : 590381833 : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1923 : : end do
1924 : 3770428 : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1925 : 3770428 : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1926 : 3770428 : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1927 : :
1928 : : case(4)
1929 : :
1930 : : scr1=0.0d0 ; sci1=0.0d0
1931 : : scr2=0.0d0 ; sci2=0.0d0
1932 : : scr3=0.0d0 ; sci3=0.0d0
1933 : : scr4=0.0d0 ; sci4=0.0d0
1934 [ + + ]: 65679398639 : do ipw=1,nincpw
1935 : 65237645501 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1936 : 65237645501 : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1937 : 65237645501 : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1938 : 65237645501 : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1939 : 65679398639 : scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
1940 : : end do
1941 : 441753138 : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1942 : 441753138 : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1943 : 441753138 : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1944 : 441753138 : scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
1945 : :
1946 : : case(5)
1947 : :
1948 : : scr1=0.0d0 ; sci1=0.0d0
1949 : : scr2=0.0d0 ; sci2=0.0d0
1950 : : scr3=0.0d0 ; sci3=0.0d0
1951 : : scr4=0.0d0 ; sci4=0.0d0
1952 : : scr5=0.0d0 ; sci5=0.0d0
1953 : : do ipw=1,nincpw
1954 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1955 : : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1956 : : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1957 : : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1958 : : scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
1959 : : scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
1960 : : end do
1961 : : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1962 : : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1963 : : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1964 : : scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
1965 : : scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
1966 : :
1967 : : case(6)
1968 : :
1969 : : scr1=0.0d0 ; sci1=0.0d0
1970 : : scr2=0.0d0 ; sci2=0.0d0
1971 : : scr3=0.0d0 ; sci3=0.0d0
1972 : : scr4=0.0d0 ; sci4=0.0d0
1973 : : scr5=0.0d0 ; sci5=0.0d0
1974 : : scr6=0.0d0 ; sci6=0.0d0
1975 : : do ipw=1,nincpw
1976 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1977 : : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1978 : : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1979 : : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1980 : : scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
1981 : : scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
1982 : : scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw)
1983 : : end do
1984 : : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1985 : : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1986 : : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1987 : : scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
1988 : : scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
1989 : : scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6
1990 : :
1991 : : case(7)
1992 : :
1993 : : scr1=0.0d0 ; sci1=0.0d0
1994 : : scr2=0.0d0 ; sci2=0.0d0
1995 : : scr3=0.0d0 ; sci3=0.0d0
1996 : : scr4=0.0d0 ; sci4=0.0d0
1997 : : scr5=0.0d0 ; sci5=0.0d0
1998 : : scr6=0.0d0 ; sci6=0.0d0
1999 : : scr7=0.0d0 ; sci7=0.0d0
2000 : : do ipw=1,nincpw
2001 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2002 : : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
2003 : : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
2004 : : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
2005 : : scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
2006 : : scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
2007 : : scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw)
2008 : : scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw)
2009 : : end do
2010 : : scalars(1,start )=scr1 ; scalars(2,start )=sci1
2011 : : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
2012 : : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
2013 : : scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
2014 : : scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
2015 : : scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6
2016 : : scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7
2017 : :
2018 : : case(8)
2019 : :
2020 : : scr1=0.0d0 ; sci1=0.0d0
2021 : : scr2=0.0d0 ; sci2=0.0d0
2022 : : scr3=0.0d0 ; sci3=0.0d0
2023 : : scr4=0.0d0 ; sci4=0.0d0
2024 : : scr5=0.0d0 ; sci5=0.0d0
2025 : : scr6=0.0d0 ; sci6=0.0d0
2026 : : scr7=0.0d0 ; sci7=0.0d0
2027 : : scr8=0.0d0 ; sci8=0.0d0
2028 : : do ipw=1,nincpw
2029 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2030 : : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
2031 : : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
2032 : : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
2033 : : scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
2034 : : scr5=scr5+ar*ffkg(start+4,ipw) ; sci5=sci5+ai*ffkg(start+4,ipw)
2035 : : scr6=scr6+ar*ffkg(start+5,ipw) ; sci6=sci6+ai*ffkg(start+5,ipw)
2036 : : scr7=scr7+ar*ffkg(start+6,ipw) ; sci7=sci7+ai*ffkg(start+6,ipw)
2037 : : scr8=scr8+ar*ffkg(start+7,ipw) ; sci8=sci8+ai*ffkg(start+7,ipw)
2038 : : end do
2039 : : scalars(1,start )=scr1 ; scalars(2,start )=sci1
2040 : : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
2041 : : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
2042 : : scalars(1,start+3)=scr4 ; scalars(2,start+3)=sci4
2043 : : scalars(1,start+4)=scr5 ; scalars(2,start+4)=sci5
2044 : : scalars(1,start+5)=scr6 ; scalars(2,start+5)=sci6
2045 : : scalars(1,start+6)=scr7 ; scalars(2,start+6)=sci7
2046 [ + + + + : 483141432 : scalars(1,start+7)=scr8 ; scalars(2,start+7)=sci8
- ]
2047 : :
2048 : : end select
2049 : :
2050 : : ! End loop on start
2051 : : end do
2052 : :
2053 : : ! End if statement for small or big nffkg
2054 : : end if
2055 : :
2056 : : ! ******* Leaving the critical part *********************************
2057 : :
2058 : : ! DEBUG
2059 : : ! write(std_out,*)' opernl4a, write scalars '
2060 : : ! do iffkg=1,nffkg
2061 : : ! write(std_out,*)iffkg,scalars(1:2,iffkg)
2062 : : ! end do
2063 : : ! ENDDEBUG
2064 : :
2065 [ + + ]: 147308012 : if(istwf_k>=2)then
2066 : : ! Impose parity of resulting scalar (this operation could be
2067 : : ! replaced by direct saving of CPU time in the preceeding section)
2068 [ + + ]: 277190368 : do iffkg=1,nffkg
2069 : 277190368 : scalars(parity(iffkg),iffkg)=0.0d0
2070 : : end do
2071 : : end if
2072 : :
2073 : 147308012 : iffkg=0 ; iffkgs=nffkge+nffkgd ; iffkgk=nffkge*2
2074 : 147308012 : iffkgs2=nffkge+nffkgs
2075 [ + + ]: 603553638 : do ilang=1,nlang
2076 : 363107250 : nproj=jproj(ilang)
2077 [ + + ]: 510415262 : if(nproj>0)then
2078 : : ! ilang2 is the number of independent tensor components
2079 : : ! for symmetric tensor of rank ilang-1
2080 : 344741790 : ilang2=(ilang*(ilang+1))/2
2081 : :
2082 : : ! Loop over projectors
2083 [ + + ]: 993824322 : do iproj=1,nproj
2084 : : ! Multiply by the k+G factors (tensors of various rank)
2085 [ + + ]: 2609501947 : do ii=1,ilang2
2086 : : ! Get the starting address for the relevant tensor
2087 : 1960419415 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
2088 : 1960419415 : iffkg=iffkg+1
2089 : : ! !$OMP CRITICAL (OPERNL4a_1)
2090 : 1960419415 : gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
2091 : 1960419415 : gxa(2,jj,ia,iproj)=gxa(2,jj,ia,iproj)+scalars(2,iffkg)
2092 : : ! !$OMP END CRITICAL (OPERNL4a_1)
2093 : : ! Now, compute gradients, if needed.
2094 [ + + + + ]: 1960419415 : if ((choice==2.or.choice==23) .and. ndgxdt==3) then
2095 [ + + ]: 101382728 : do mu=1,3
2096 : 76037046 : jffkg=nffkge+(iffkg-1)*3+mu
2097 : : ! Pay attention to the use of reals and imaginary parts here ...
2098 : : ! !$OMP CRITICAL (OPERNL4a_2)
2099 : 76037046 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2100 : 101382728 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
2101 : : ! !$OMP END CRITICAL (OPERNL4a_2)
2102 : : end do
2103 : : end if
2104 [ + + + + ]: 1960419415 : if (choice==2 .and. ndgxdt==1) then
2105 : 144190839 : jffkg=nffkge+iffkg
2106 : : ! Pay attention to the use of reals and imaginary parts here ...
2107 : : ! !$OMP CRITICAL (OPERNL4a_3)
2108 : 144190839 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2109 : 144190839 : dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+two_pi*scalars(1,jffkg)
2110 : : ! !$OMP END CRITICAL (OPERNL4a_3)
2111 : : end if
2112 [ + + ]: 2609501947 : if (choice==4) then
2113 [ + + ]: 23749800 : do mu=1,3
2114 : 17812350 : jffkg=nffkge+(iffkg-1)*9+mu
2115 : : ! Pay attention to the use of reals and imaginary parts here ...
2116 : : ! !$OMP CRITICAL (OPERNL4a_4)
2117 : 17812350 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2118 : 23749800 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)+two_pi*scalars(1,jffkg)
2119 : : ! !$OMP END CRITICAL (OPERNL4a_4)
2120 : : end do
2121 [ + + ]: 41562150 : do mu=4,9
2122 : 35624700 : jffkg=nffkge+(iffkg-1)*9+mu
2123 : : ! Pay attention to the use of reals and imaginary parts here ...
2124 : : ! Also, note the multiplication by (2 pi)**2
2125 : : ! !$OMP CRITICAL (OPERNL4a_5)
2126 : 35624700 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi2*scalars(1,jffkg)
2127 : 41562150 : dgxdt(2,mu,jj,ia,iproj)=dgxdt(2,mu,jj,ia,iproj)-two_pi2*scalars(2,jffkg)
2128 : : ! !$OMP END CRITICAL (OPERNL4a_5)
2129 : : end do
2130 : : end if
2131 : : ! End loop on ii=1,ilang2
2132 : : end do
2133 : :
2134 : : if ((choice==3.or.choice==23) .or. choice==6) then
2135 : : ! Compute additional tensors related to strain gradients
2136 : : ! ilang3 is number of unique tensor components of rank ilang+1
2137 : 10595410 : ilang3=((ilang+2)*(ilang+3))/2
2138 : 10595410 : jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
2139 : : ! Compute strain gradient tensor components
2140 [ + + ]: 119435029 : do ii=1,ilang3
2141 : : ! Note that iffkgs is also used by ddk and 2nd derivative parts
2142 : 108839619 : iffkgs=iffkgs+1
2143 : 108839619 : jj=ii+jjs
2144 : : ! !$OMP CRITICAL (OPERNL4a_6)
2145 : 108839619 : dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
2146 : 119435029 : dgxds(2,jj-4,ia,iproj)=dgxds(2,jj-4,ia,iproj)+scalars(2,iffkgs)
2147 : : ! !$OMP END CRITICAL (OPERNL4a_6)
2148 : : end do
2149 : : end if
2150 : :
2151 [ + + ]: 649082532 : if (choice==6) then
2152 : : ! Compute additional tensors related to strain 2nd derivatives
2153 : : ! and internal strain derivatives
2154 : : ! ilang6 is number of unique tensor components of rank ilang+3
2155 : 179208 : ilang6=((ilang+4)*(ilang+5))/2
2156 : 179208 : jjs=((ilang+3)*(ilang+4)*(ilang+5))/6
2157 : : ! Compute strain gradient tensor components
2158 [ + + ]: 3973648 : do ii=1,ilang6
2159 : : ! Note that iffkgs is also used by ddk part
2160 : 3794440 : iffkgs2=iffkgs2+1
2161 : 3794440 : jj=ii+jjs
2162 : : ! !$OMP CRITICAL (OPERNL4a_7)
2163 : 3794440 : d2gxds2(1,jj-20,ia,iproj)=d2gxds2(1,jj-20,ia,iproj)+scalars(1,iffkgs2)
2164 : 3973648 : d2gxds2(2,jj-20,ia,iproj)=d2gxds2(2,jj-20,ia,iproj)+scalars(2,iffkgs2)
2165 : : ! !$OMP END CRITICAL (OPERNL4a_7)
2166 : : end do
2167 : :
2168 : : ! ilang4 is number of unique tensor components of rank ilang
2169 : 179208 : ilang4=((ilang+1)*(ilang+2))/2
2170 : 179208 : jjs=((ilang)*(ilang+1)*(ilang+2))/6
2171 : : ! Compute internal strain gradient tensor components
2172 [ + + ]: 1298728 : do ii=1,ilang4
2173 : 1119520 : iffkgs2=iffkgs2+1
2174 : 1119520 : jj=ii+jjs
2175 : : ! !$OMP CRITICAL (OPERNL4a_8)
2176 : : ! Pay attention to the use of reals and imaginary parts here ...
2177 : 1119520 : dgxdis(1,jj-1,ia,iproj)=dgxdis(1,jj-1,ia,iproj)-two_pi*scalars(2,iffkgs2)
2178 : 1298728 : dgxdis(2,jj-1,ia,iproj)=dgxdis(2,jj-1,ia,iproj)+two_pi*scalars(1,iffkgs2)
2179 : : ! !$OMP END CRITICAL (OPERNL4a_8)
2180 : : end do
2181 : :
2182 : : ! ilang5 is number of unique tensor components of rank ilang+2
2183 : 179208 : ilang5=((ilang+3)*(ilang+4))/2
2184 : 179208 : jjs=((ilang+2)*(ilang+3)*(ilang+4))/6
2185 : : ! Compute internal strain gradient tensor components
2186 [ + + ]: 2902800 : do ii=1,ilang5
2187 : 2723592 : iffkgs2=iffkgs2+1
2188 : 2723592 : jj=ii+jjs
2189 : : ! !$OMP CRITICAL (OPERNL4a_9)
2190 : : ! Pay attention to the use of reals and imaginary parts here ...
2191 : 2723592 : d2gxdis(1,jj-10,ia,iproj)=d2gxdis(1,jj-10,ia,iproj)-two_pi*scalars(2,iffkgs2)
2192 : 2902800 : d2gxdis(2,jj-10,ia,iproj)=d2gxdis(2,jj-10,ia,iproj)+two_pi*scalars(1,iffkgs2)
2193 : : ! !$OMP END CRITICAL (OPERNL4a_9)
2194 : : end do
2195 : : end if ! choice==6
2196 : :
2197 [ + + ]: 993824322 : if (choice==5) then
2198 : : ! Compute additional tensors related to ddk with ffnl(:,2,...)
2199 : 7800131 : ilangx=(ilang*(ilang+1))/2
2200 : 7800131 : jjs=((ilang-1)*ilang*(ilang+1))/6
2201 [ + + ]: 31880060 : do ii=1,ilangx
2202 : : ! Note that iffkgs is also used by strain part
2203 : 24079929 : iffkgs=iffkgs+1
2204 : 24079929 : jj=ii+jjs
2205 : : ! !$OMP CRITICAL (OPERNL4a_10)
2206 : 24079929 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
2207 : 31880060 : dgxdt(2,1,jj,ia,iproj)=dgxdt(2,1,jj,ia,iproj)+scalars(2,iffkgs)
2208 : : ! !$OMP END CRITICAL (OPERNL4a_10)
2209 : : end do
2210 : : ! Compute additional tensors related to ddk with ffnl(:,1,...)
2211 [ + + ]: 7800131 : if(ilang>=2)then
2212 : 4980659 : ilangx=((ilang-1)*ilang)/2
2213 : 4980659 : jjs=((ilang-2)*(ilang-1)*ilang)/6
2214 [ + + ]: 14174778 : do ii=1,ilangx
2215 : 9194119 : iffkgk=iffkgk+1
2216 : 9194119 : jj=ii+jjs
2217 : : ! !$OMP CRITICAL (OPERNL4a_11)
2218 : 9194119 : dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
2219 : 14174778 : dgxdt(2,2,jj,ia,iproj)=dgxdt(2,2,jj,ia,iproj)+scalars(2,iffkgk)
2220 : : ! !$OMP END CRITICAL (OPERNL4a_11)
2221 : : end do
2222 : : end if
2223 : : end if
2224 : :
2225 : : ! End projector loop
2226 : : end do
2227 : :
2228 : : ! End condition of non-zero projectors
2229 : : end if
2230 : :
2231 : : ! End angular momentum loop
2232 : : end do
2233 : :
2234 : : ! End loop on atoms
2235 : : end do
2236 : :
2237 : : ! End loop on blocks of planewaves
2238 : : end do
2239 : : !!$OMP END DO
2240 : :
2241 [ - + ]: 48511395 : ABI_FREE(ffkg)
2242 [ - + ]: 48511395 : ABI_FREE(kpgx)
2243 [ - + ]: 48511395 : ABI_FREE(parity)
2244 [ - + ]: 48511395 : ABI_FREE(scalars)
2245 [ - + ]: 48511395 : ABI_FREE(teffv)
2246 : : !!$OMP END PARALLEL
2247 : :
2248 : : !DEBUG
2249 : : !write(std_out,*)' opernl4a : exit'
2250 : : !ENDDEBUG
2251 : :
2252 : 48511395 : end subroutine opernl4a
2253 : : !!***
2254 : :
2255 : : !!****f* ABINIT/opernl4b
2256 : : !! NAME
2257 : : !! opernl4b
2258 : : !!
2259 : : !! FUNCTION
2260 : : !! Operate with the non-local part of the hamiltonian,
2261 : : !! from projected quantities to reciprocal space
2262 : : !!
2263 : : !! INPUTS
2264 : : !! if(choice==2 .or choice==4 .or. choice==5)
2265 : : !! dgxdt(2,ndgxdt,mlang3,mincat,mproj)= selected gradients of gxa wrt coords
2266 : : !! or with respect to ddk
2267 : : !! if(choice==3)
2268 : : !! dgxds((2,mlang4,mincat,mproj) = gradients of projected scalars wrt strains
2269 : : !! ------ Taken away in beautification because unused MS -------
2270 : : !! d2gxds2((2,mlang6,mincat,mproj) dummy argument here, not used
2271 : : !! -------------------------------------------------------------
2272 : : !! ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
2273 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
2274 : : !! gxa(2,mlang3,mincat,mproj)= projected scalars
2275 : : !! ia3=gives the number of the first atom in the subset presently treated
2276 : : !! idir=direction of the perturbation (needed if choice==2 or 5, and ndgxdt=1)
2277 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln
2278 : : !! ispinor=1 or 2, gives the spinorial component of ffnl to be used
2279 : : !! ------ Taken away in beautification because unused MS -------
2280 : : !! istwf_k=option parameter that describes the storage of wfs
2281 : : !! -------------------------------------------------------------
2282 : : !! itypat = type of atom, needed for ffnl
2283 : : !! jproj(nlang)=number of projectors for each angular momentum
2284 : : !! kg_k(3,npw)=integer coords of planewaves in basis sphere
2285 : : !! kpg_k(npw,npkg)= (k+G) components and related data
2286 : : !! kpt(3)=real components of k point in terms of recip. translations
2287 : : !! lmnmax=max. number of (l,n) components over all type of psps
2288 : : !! matblk=dimension of the array ph3d
2289 : : !! mincat= maximum increment of atoms
2290 : : !! mlang3 = one of the dimensions of the array gxa
2291 : : !! mlang4 = dimension for dgxds
2292 : : !! ------ Taken away in beautification because unused MS -------
2293 : : !! mlang6 = dimension for d2gxds2
2294 : : !! -------------------------------------------------------------
2295 : : !! mproj=maximum dimension for number of projection operators for each
2296 : : !! angular momentum for nonlocal pseudopotential
2297 : : !! ndgxdt=second dimension of dgxdt
2298 : : !! nincat = number of atoms in the subset here treated
2299 : : !! nkpg=second size of array kpg_k
2300 : : !! nffnl=third dimension of ffnl
2301 : : !! nlang = number of angular momenta to be treated = 1 + highest ang. mom.
2302 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
2303 : : !! npw = number of plane waves in reciprocal space
2304 : : !! ntypat = number of type of atoms, dimension needed for ffnl
2305 : : !! gxa(2,mlang3,nincat,mproj)=modified projected scalars;
2306 : : !! NOTE that metric contractions have already been performed
2307 : : !! ph3d(2,npw,matblk)=three-dimensional phase factors
2308 : : !!
2309 : : !! OUTPUT
2310 : : !! vect(2*npw)=final vector in reciprocal space <G|V_nonlocal|vect_start>.
2311 : : !!
2312 : : !! NOTES
2313 : : !! Operate with the non-local part of the hamiltonian for one type of
2314 : : !! atom, and within this given type of atom, for a subset of
2315 : : !! at most nincat atoms.
2316 : : !!
2317 : : !! This routine basically replaces getgla (gxa here is the former gla),
2318 : : !! except for the calculation of <G|dVnl/dk|C> or strain gradients.
2319 : : !!
2320 : : !! Present version decomposed according to iffkg
2321 : : !! opernl4a.f is from reciprocal space to projected quantities.
2322 : : !!
2323 : : !! SOURCE
2324 : :
2325 : 47822681 : subroutine opernl4b(choice,dgxds,dgxdt,ffnl,gmet,gxa,&
2326 : 47822681 : & ia3,idir,indlmn,ispinor,itypat,jproj,kg_k,kpg_k,kpt,&
2327 : : & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat,&
2328 [ - + ]: 47822681 : & nkpg,nlang,nloalg,npw,ntypat,ph3d,vect)
2329 : :
2330 : : !Arguments ------------------------------------
2331 : : !scalars
2332 : : integer,intent(in) :: choice,ia3,idir,ispinor,itypat,lmnmax,matblk !,istwf_k
2333 : : integer,intent(in) :: mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat !,mlang6
2334 : : integer,intent(in) :: nkpg,nlang,npw,ntypat
2335 : : !arrays
2336 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),jproj(nlang),kg_k(3,npw)
2337 : : integer,intent(in) :: nloalg(3)
2338 : : !real(dp),intent(in) :: d2gxds2(2,mlang6,mincat,mproj)
2339 : : real(dp),intent(in) :: dgxds(2,mlang4,mincat,mproj)
2340 : : real(dp),intent(in) :: dgxdt(2,ndgxdt,mlang3,mincat,mproj)
2341 : : real(dp),intent(in) :: ffnl(1,npw,nffnl,lmnmax,ntypat),gmet(3,3)
2342 : : real(dp),intent(in) :: gxa(2,mlang3,mincat,mproj),kpg_k(npw,nkpg),kpt(3)
2343 : : real(dp),intent(in) :: ph3d(2,npw,matblk)
2344 : : real(dp),intent(inout) :: vect(:,:) !vz_i
2345 : :
2346 : : !Local variables-------------------------------
2347 : : !scalars
2348 : : integer :: chunk,ia,iaph3d,iffkg,iffkgk,iffkgs,ig,ii,ilang,ilang2,ilang3
2349 : : integer :: iproj,ipw,ipw1,ipw2,jffkg,jj,jump,mblkpw,nffkg
2350 : : integer :: nffkgd,nffkge,nffkgk,nffkgs,nincpw,nproj,ntens,start
2351 : : real(dp) :: ai,ar,sci1,sci2,sci3,sci4,sci5,sci6,sci7,sci8
2352 : : real(dp) :: scr1,scr2,scr3,scr4,scr5,scr6,scr7,scr8
2353 : : character(len=500) :: message
2354 : : !arrays
2355 : 47822681 : integer,allocatable :: parity(:)
2356 : 47822681 : real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
2357 : :
2358 : : ! *************************************************************************
2359 : :
2360 : : !mblkpw sets the size of blocks of planewaves
2361 : 47822681 : mblkpw=NLO_MBLKPW
2362 : :
2363 : : !jump governs, in fine, the use of registers in the most cpu
2364 : : !time consuming part of the routine. Until now, jump=8 is the maximal value.
2365 : : !The optimal value will be machine-dependent !
2366 : 47822681 : jump=4
2367 : :
2368 : : !Initialisation before blocking on the plane waves
2369 : :
2370 : : !Set up dimension of kpgx and allocate
2371 : : !ntens sets the maximum number of independent tensor components
2372 : : !over all allowed angular momenta; need 20 for spdf for tensors
2373 : : !up to rank 3; to handle stress tensor, need up to rank 5
2374 : 47822681 : ntens=1
2375 [ + + + + : 47822681 : if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5)ntens=4
+ - + + ]
2376 [ + + + + ]: 47822681 : if(nlang>=3 .or. choice==3)ntens=10
2377 [ + + + + : 47822681 : if(nlang>=4 .or. (choice==3 .and. nlang>=2) )ntens=20
+ + ]
2378 [ + + + + ]: 47822681 : if(choice==3 .and. nlang>=3)ntens=35
2379 [ + + - + ]: 47822681 : if(choice==3 .and. nlang==4)ntens=56
2380 : :
2381 : : !Set up second dimension of ffkg array, and allocate
2382 : 47822681 : nffkg=0; nffkge=0; nffkgd=0; nffkgk=0; nffkgs=0
2383 [ + + ]: 164692855 : do ilang=1,nlang
2384 : : ! Get the number of projectors for that angular momentum
2385 : 116870174 : nproj=jproj(ilang)
2386 : : ! If there is a non-local part, accumulate the number of vectors needed
2387 [ + + ]: 164692855 : if(nproj>0)then
2388 : 110414140 : ilang2=(ilang*(ilang+1))/2
2389 : 110414140 : nffkge=nffkge+nproj*ilang2
2390 [ + + ]: 110414140 : if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
2391 [ + + ]: 110414140 : if(choice==2 .or. choice==4)nffkgd=nffkgd+ndgxdt*nproj*ilang2
2392 [ + + ]: 110414140 : if(choice==3)then
2393 : 2353575 : ilang3=((ilang+2)*(ilang+3))/2
2394 : 2353575 : nffkgs=nffkgs+nproj*ilang3
2395 : : end if
2396 : : end if
2397 : : end do
2398 : 47822681 : nffkg=nffkge+nffkgd+nffkgs+nffkgk
2399 : :
2400 : : !!$OMP PARALLEL DEFAULT(PRIVATE) &
2401 : : !!$OMP SHARED(choice,dgxds,dgxdt,ffnl,gmet,gxa,ia3,idir,indlmn,ispinor) &
2402 : : !!$OMP SHARED(itypat,jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang3,mlang4,mproj) &
2403 : : !!$OMP SHARED(ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,ntypat,ph3d,vect) &
2404 : : !!$OMP SHARED(jump,nffkgd,nffkgk,nffkgs,mblkpw,nffkg,nffkge,ntens)
2405 : :
2406 [ - + + - : 191290724 : ABI_MALLOC(ffkg,(nffkg,mblkpw))
+ - + - +
- - + -
+ ]
2407 [ - + + - : 143468043 : ABI_MALLOC(parity,(nffkg))
+ - + - -
+ - + ]
2408 [ - + + - : 143468043 : ABI_MALLOC(kpgx,(mblkpw,ntens))
+ - + - -
+ - + ]
2409 [ - + + - : 143468043 : ABI_MALLOC(scalars,(2,nffkg))
+ - + - -
+ - + ]
2410 [ - + ]: 47822681 : ABI_MALLOC(teffv,(2,mblkpw))
2411 : :
2412 : : !Loop on subsets of plane waves (blocking)
2413 : : !!$OMP DO
2414 [ + - + + ]: 139559355 : do ipw1=1,npw,mblkpw
2415 : :
2416 : 91736674 : ipw2=min(npw,ipw1+mblkpw-1)
2417 : 91736674 : nincpw=ipw2-ipw1+1
2418 : :
2419 : : ! Initialize kpgx array related to tensors defined below
2420 : : call dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
2421 : : & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
2422 : 91736674 : & npw,ntens,ntypat,parity)
2423 : :
2424 [ + - ]: 139559355 : if (choice==1 .or. choice==2 .or. choice==3 .or. choice==5) then
2425 : : ! Application of non-local part from projected scalars
2426 : : ! back to reciprocal space ...
2427 : : ! [this section merely computes terms which add to <G|Vnl|C>;
2428 : : ! nothing here is needed when various gradients are being computed]
2429 : :
2430 : : ! Loop on atoms
2431 [ + + ]: 234777441 : do ia=1,nincat
2432 : :
2433 : : ! Compute the shift eventually needed to get the phases in ph3d
2434 : 143040767 : iaph3d=ia
2435 [ + + ]: 143040767 : if(nloalg(2)>0)iaph3d=ia+ia3-1
2436 : :
2437 : : ! Transfer gxa (and eventually dgxdt) in scalars with different indexing
2438 : 143040767 : iffkg=0
2439 : 143040767 : iffkgk=nffkge*2
2440 : 143040767 : iffkgs=nffkge
2441 [ + + ]: 495014985 : do ilang=1,nlang
2442 : 351974218 : nproj=jproj(ilang)
2443 [ + + ]: 495014985 : if (nproj>0) then
2444 : 333866529 : ilang2=(ilang*(ilang+1))/2
2445 : 333866529 : ilang3=((ilang+2)*(ilang+3))/2
2446 [ + + ]: 962053614 : do iproj=1,nproj
2447 [ + + ]: 2522062762 : do ii=1,ilang2
2448 : 1893875677 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
2449 : 1893875677 : iffkg=iffkg+1
2450 [ + + ]: 2522062762 : if(choice==1 .or. choice==3)then
2451 : 1725352589 : scalars(1,iffkg)=gxa(1,jj,ia,iproj)
2452 : 1725352589 : scalars(2,iffkg)=gxa(2,jj,ia,iproj)
2453 [ + + + - ]: 168523088 : else if (choice==2 .and. ndgxdt==1) then
2454 : 144443159 : jffkg=nffkge+iffkg
2455 : : ! Pay attention to the use of reals and imaginary parts here ...
2456 : : ! Also, the gxa and dgxdt arrays are switched, in order
2457 : : ! to give the correct combination when multiplying ffkg,
2458 : : ! see Eq.(53) of PRB55,10337(1997) [[cite:Gonze1997]]
2459 : 144443159 : scalars(1,jffkg)= two_pi*gxa(2,jj,ia,iproj)
2460 : 144443159 : scalars(2,jffkg)=-two_pi*gxa(1,jj,ia,iproj)
2461 : 144443159 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
2462 : 144443159 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
2463 [ + - ]: 24079929 : else if (choice==5) then
2464 : 24079929 : jffkg=nffkge+iffkg
2465 : : ! The gxa and dgxdt arrays are switched, in order
2466 : : ! to give the correct combination when multiplying ffkg,
2467 : 24079929 : scalars(1,jffkg)= gxa(1,jj,ia,iproj)
2468 : 24079929 : scalars(2,jffkg)= gxa(2,jj,ia,iproj)
2469 : 24079929 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
2470 : 24079929 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
2471 : : end if
2472 : : end do
2473 [ + + ]: 628187085 : if(choice==3) then
2474 [ + + ]: 58947366 : do ii=1,ilang3
2475 : 53697342 : iffkgs=iffkgs+1
2476 : 53697342 : jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
2477 : 53697342 : scalars(1,iffkgs)=dgxds(1,jj-4,ia,iproj)
2478 : 58947366 : scalars(2,iffkgs)=dgxds(2,jj-4,ia,iproj)
2479 : : end do
2480 : : end if
2481 [ + + + + ]: 962053614 : if(ilang>=2 .and. choice==5)then
2482 [ + + ]: 14174778 : do ii=1,((ilang-1)*ilang)/2
2483 : 9194119 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
2484 : 9194119 : iffkgk=iffkgk+1
2485 : 9194119 : scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
2486 : 14174778 : scalars(2,iffkgk)= dgxdt(2,2,jj,ia,iproj)
2487 : : end do
2488 : : end if
2489 : : end do
2490 : : end if
2491 : : end do
2492 : :
2493 : : ! DEBUG
2494 : : ! write(std_out,*)' opernl4b, write scalars '
2495 : : ! do iffkg=1,nffkg
2496 : : ! write(std_out,*)iffkg,scalars(1:2,iffkg)
2497 : : ! end do
2498 : : ! ENDDEBUG
2499 : :
2500 : : ! ******* Entering the second critical part ****************************
2501 : :
2502 : : ! First, treat small nffkg; send treat the loop needed for big nffkg;
2503 : : ! finally treat the end of the loop needed for big nffkg
2504 : :
2505 : : ! For the time being, the maximal jump allowed is 8.
2506 : :
2507 : : ! 1) Here, treat small nffkg
2508 [ + + ]: 234777441 : if(nffkg<=jump)then
2509 : :
2510 : 13057071 : select case(nffkg)
2511 : :
2512 : : case(1)
2513 : :
2514 : 13057071 : scr1=scalars(1,1) ; sci1=scalars(2,1)
2515 : 13057071 : ig=ipw1
2516 [ + + ]: 1791610038 : do ipw=1,nincpw
2517 : 1778552967 : ar=ffkg(1,ipw)*scr1 ; ai=ffkg(1,ipw)*sci1
2518 : 1778552967 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2519 : 1778552967 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2520 : 1791610038 : ig=ig+1
2521 : : end do
2522 : :
2523 : : case(2)
2524 : :
2525 : 2362732 : scr1=scalars(1,1) ; sci1=scalars(2,1)
2526 : 2362732 : scr2=scalars(1,2) ; sci2=scalars(2,2)
2527 : 2362732 : ig=ipw1
2528 [ + + ]: 320123169 : do ipw=1,nincpw
2529 : 317760437 : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2530 : 317760437 : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2531 : 317760437 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2532 : 317760437 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2533 : 320123169 : ig=ig+1
2534 : : end do
2535 : :
2536 : : case(3)
2537 : :
2538 : 6856 : scr1=scalars(1,1) ; sci1=scalars(2,1)
2539 : 6856 : scr2=scalars(1,2) ; sci2=scalars(2,2)
2540 : 6856 : scr3=scalars(1,3) ; sci3=scalars(2,3)
2541 : 6856 : ig=ipw1
2542 [ + + ]: 1359388 : do ipw=1,nincpw
2543 : 1352532 : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2544 : 1352532 : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2545 : 1352532 : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2546 : 1352532 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2547 : 1352532 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2548 : 1359388 : ig=ig+1
2549 : : end do
2550 : :
2551 : : case(4)
2552 : :
2553 : 12629261 : scr1=scalars(1,1) ; sci1=scalars(2,1)
2554 : 12629261 : scr2=scalars(1,2) ; sci2=scalars(2,2)
2555 : 12629261 : scr3=scalars(1,3) ; sci3=scalars(2,3)
2556 : 12629261 : scr4=scalars(1,4) ; sci4=scalars(2,4)
2557 : 12629261 : ig=ipw1
2558 [ + + ]: 1566714450 : do ipw=1,nincpw
2559 : 1554085189 : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2560 : 1554085189 : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2561 : 1554085189 : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2562 : 1554085189 : ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2563 : 1554085189 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2564 : 1554085189 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2565 : 1566714450 : ig=ig+1
2566 : : end do
2567 : :
2568 : : case(5)
2569 : :
2570 : : scr1=scalars(1,1) ; sci1=scalars(2,1)
2571 : : scr2=scalars(1,2) ; sci2=scalars(2,2)
2572 : : scr3=scalars(1,3) ; sci3=scalars(2,3)
2573 : : scr4=scalars(1,4) ; sci4=scalars(2,4)
2574 : : scr5=scalars(1,5) ; sci5=scalars(2,5)
2575 : : ig=ipw1
2576 : : do ipw=1,nincpw
2577 : : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2578 : : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2579 : : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2580 : : ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2581 : : ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2582 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2583 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2584 : : ig=ig+1
2585 : : end do
2586 : :
2587 : : case(6)
2588 : :
2589 : : scr1=scalars(1,1) ; sci1=scalars(2,1)
2590 : : scr2=scalars(1,2) ; sci2=scalars(2,2)
2591 : : scr3=scalars(1,3) ; sci3=scalars(2,3)
2592 : : scr4=scalars(1,4) ; sci4=scalars(2,4)
2593 : : scr5=scalars(1,5) ; sci5=scalars(2,5)
2594 : : scr6=scalars(1,6) ; sci6=scalars(2,6)
2595 : : ig=ipw1
2596 : : do ipw=1,nincpw
2597 : : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2598 : : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2599 : : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2600 : : ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2601 : : ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2602 : : ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
2603 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2604 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2605 : : ig=ig+1
2606 : : end do
2607 : :
2608 : : case(7)
2609 : :
2610 : : scr1=scalars(1,1) ; sci1=scalars(2,1)
2611 : : scr2=scalars(1,2) ; sci2=scalars(2,2)
2612 : : scr3=scalars(1,3) ; sci3=scalars(2,3)
2613 : : scr4=scalars(1,4) ; sci4=scalars(2,4)
2614 : : scr5=scalars(1,5) ; sci5=scalars(2,5)
2615 : : scr6=scalars(1,6) ; sci6=scalars(2,6)
2616 : : scr7=scalars(1,7) ; sci7=scalars(2,7)
2617 : : ig=ipw1
2618 : : do ipw=1,nincpw
2619 : : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2620 : : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2621 : : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2622 : : ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2623 : : ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2624 : : ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
2625 : : ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7
2626 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2627 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2628 : : ig=ig+1
2629 : : end do
2630 : :
2631 : : case(8)
2632 : :
2633 : : scr1=scalars(1,1) ; sci1=scalars(2,1)
2634 : : scr2=scalars(1,2) ; sci2=scalars(2,2)
2635 : : scr3=scalars(1,3) ; sci3=scalars(2,3)
2636 : : scr4=scalars(1,4) ; sci4=scalars(2,4)
2637 : : scr5=scalars(1,5) ; sci5=scalars(2,5)
2638 : : scr6=scalars(1,6) ; sci6=scalars(2,6)
2639 : : scr7=scalars(1,7) ; sci7=scalars(2,7)
2640 : : scr8=scalars(1,8) ; sci8=scalars(2,8)
2641 : : ig=ipw1
2642 [ + + + + : 28055920 : do ipw=1,nincpw
- ]
2643 : : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2644 : : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2645 : : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2646 : : ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2647 : : ar=ar+ffkg(5,ipw)*scr5 ; ai=ai+ffkg(5,ipw)*sci5
2648 : : ar=ar+ffkg(6,ipw)*scr6 ; ai=ai+ffkg(6,ipw)*sci6
2649 : : ar=ar+ffkg(7,ipw)*scr7 ; ai=ai+ffkg(7,ipw)*sci7
2650 : : ar=ar+ffkg(8,ipw)*scr8 ; ai=ai+ffkg(8,ipw)*sci8
2651 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2652 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2653 : : ig=ig+1
2654 : : end do
2655 : :
2656 : :
2657 : : end select
2658 : :
2659 : : else
2660 : :
2661 [ + + ]: 18037739031 : do ipw=1,nincpw
2662 : 18037739031 : teffv(1,ipw)=0.0d0 ; teffv(2,ipw)=0.0d0
2663 : : end do
2664 : :
2665 : : ! 2) Here treart the loop for big nffkg
2666 [ + - + + ]: 537254819 : do start=1,nffkg-jump,jump
2667 : 422269972 : chunk=min(jump,nffkg-jump-start+1)
2668 : :
2669 : 114984847 : select case(chunk)
2670 : :
2671 : : case(1)
2672 : :
2673 : 14964285 : scr1=scalars(1,start) ; sci1=scalars(2,start)
2674 [ + + ]: 2133681487 : do ipw=1,nincpw
2675 : 2118717202 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2676 : 2118717202 : ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1
2677 : 2133681487 : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2678 : : end do
2679 : :
2680 : : case(2)
2681 : :
2682 : 21758230 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2683 : 21758230 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2684 [ + + ]: 3520018343 : do ipw=1,nincpw
2685 : 3498260113 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2686 : 3498260113 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2687 : 3498260113 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2688 : 3520018343 : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2689 : : end do
2690 : :
2691 : : case(3)
2692 : :
2693 : 3639882 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2694 : 3639882 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2695 : 3639882 : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2696 [ + + ]: 570030789 : do ipw=1,nincpw
2697 : 566390907 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2698 : 566390907 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2699 : 566390907 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2700 : 566390907 : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2701 : 570030789 : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2702 : : end do
2703 : :
2704 : : case(4)
2705 : :
2706 : 381907575 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2707 : 381907575 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2708 : 381907575 : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2709 : 381907575 : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2710 [ + + ]: 56798949474 : do ipw=1,nincpw
2711 : 56417041899 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2712 : 56417041899 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2713 : 56417041899 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2714 : 56417041899 : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2715 : 56417041899 : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2716 : 56798949474 : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2717 : : end do
2718 : :
2719 : : case(5)
2720 : :
2721 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2722 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2723 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2724 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2725 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2726 : : do ipw=1,nincpw
2727 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2728 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2729 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2730 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2731 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2732 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2733 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2734 : : end do
2735 : :
2736 : : case(6)
2737 : :
2738 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2739 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2740 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2741 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2742 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2743 : : scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2744 : : do ipw=1,nincpw
2745 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2746 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2747 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2748 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2749 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2750 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2751 : : ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2752 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2753 : : end do
2754 : :
2755 : : case(7)
2756 : :
2757 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2758 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2759 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2760 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2761 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2762 : : scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2763 : : scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2764 : : do ipw=1,nincpw
2765 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2766 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2767 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2768 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2769 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2770 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2771 : : ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2772 : : ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2773 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2774 : : end do
2775 : :
2776 : : case(8)
2777 : :
2778 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2779 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2780 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2781 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2782 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2783 : : scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2784 : : scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2785 : : scr8=scalars(1,start+7) ; sci8=scalars(2,start+7)
2786 [ + + + + : 422269972 : do ipw=1,nincpw
- ]
2787 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2788 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2789 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2790 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2791 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2792 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2793 : : ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2794 : : ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2795 : : ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8
2796 : : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
2797 : : end do
2798 : :
2799 : : end select
2800 : :
2801 : : end do
2802 : :
2803 : : ! 3) Treat the end of the loops
2804 : :
2805 : 114984847 : start=nffkg-jump+1
2806 : :
2807 : : select case(jump)
2808 : :
2809 : : case(1)
2810 : :
2811 : : scr1=scalars(1,start) ; sci1=scalars(2,start)
2812 : : ig=ipw1
2813 : : do ipw=1,nincpw
2814 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2815 : : ar=ar+ffkg(start,ipw)*scr1 ; ai=ai+ffkg(start,ipw)*sci1
2816 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2817 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2818 : : ig=ig+1
2819 : : end do
2820 : :
2821 : : case(2)
2822 : :
2823 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2824 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2825 : : ig=ipw1
2826 : : do ipw=1,nincpw
2827 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2828 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2829 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2830 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2831 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2832 : : ig=ig+1
2833 : : end do
2834 : :
2835 : : case(3)
2836 : :
2837 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2838 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2839 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2840 : : ig=ipw1
2841 : : do ipw=1,nincpw
2842 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2843 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2844 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2845 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2846 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2847 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2848 : : ig=ig+1
2849 : : end do
2850 : :
2851 : : case(4)
2852 : :
2853 : 114984847 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2854 : 114984847 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2855 : 114984847 : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2856 : 114984847 : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2857 : 114984847 : ig=ipw1
2858 [ + + ]: 18037739031 : do ipw=1,nincpw
2859 : 17922754184 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2860 : 17922754184 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2861 : 17922754184 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2862 : 17922754184 : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2863 : 17922754184 : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2864 : 17922754184 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2865 : 17922754184 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2866 : 18037739031 : ig=ig+1
2867 : : end do
2868 : :
2869 : : case(5)
2870 : :
2871 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2872 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2873 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2874 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2875 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2876 : : ig=ipw1
2877 : : do ipw=1,nincpw
2878 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2879 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2880 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2881 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2882 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2883 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2884 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2885 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2886 : : ig=ig+1
2887 : : end do
2888 : :
2889 : : case(6)
2890 : :
2891 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2892 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2893 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2894 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2895 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2896 : : scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2897 : : ig=ipw1
2898 : : do ipw=1,nincpw
2899 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2900 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2901 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2902 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2903 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2904 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2905 : : ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2906 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2907 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2908 : : ig=ig+1
2909 : : end do
2910 : :
2911 : : case(7)
2912 : :
2913 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2914 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2915 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2916 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2917 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2918 : : scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2919 : : scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2920 : : ig=ipw1
2921 : : do ipw=1,nincpw
2922 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2923 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2924 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2925 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2926 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2927 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2928 : : ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2929 : : ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2930 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2931 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2932 : : ig=ig+1
2933 : : end do
2934 : :
2935 : : case(8)
2936 : :
2937 : : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2938 : : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2939 : : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2940 : : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2941 : : scr5=scalars(1,start+4) ; sci5=scalars(2,start+4)
2942 : : scr6=scalars(1,start+5) ; sci6=scalars(2,start+5)
2943 : : scr7=scalars(1,start+6) ; sci7=scalars(2,start+6)
2944 : : scr8=scalars(1,start+7) ; sci8=scalars(2,start+7)
2945 : : ig=ipw1
2946 : : do ipw=1,nincpw
2947 : : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2948 : : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2949 : : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2950 : : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2951 : : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2952 : : ar=ar+ffkg(start+4,ipw)*scr5 ; ai=ai+ffkg(start+4,ipw)*sci5
2953 : : ar=ar+ffkg(start+5,ipw)*scr6 ; ai=ai+ffkg(start+5,ipw)*sci6
2954 : : ar=ar+ffkg(start+6,ipw)*scr7 ; ai=ai+ffkg(start+6,ipw)*sci7
2955 : : ar=ar+ffkg(start+7,ipw)*scr8 ; ai=ai+ffkg(start+7,ipw)*sci8
2956 : : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2957 : : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2958 : : ig=ig+1
2959 : : end do
2960 : :
2961 : : end select
2962 : :
2963 : : ! End if statement for small or big nffkg
2964 : : end if
2965 : :
2966 : : ! ******* Leaving the critical part *********************************
2967 : :
2968 : : ! End loop on atoms
2969 : : end do
2970 : :
2971 : : ! End choice==1 or choice==2 or choice==3
2972 : : else
2973 : : ! Problem: choice does not make sense
2974 : 0 : write(message,'(a,i0,a)' )' Input choice=',choice,' not allowed. '
2975 : 0 : ABI_BUG(message)
2976 : : end if
2977 : :
2978 : : ! End loop on blocks of planewaves
2979 : : end do
2980 : : !!$OMP END DO
2981 : :
2982 [ - + ]: 47822681 : ABI_FREE(ffkg)
2983 [ - + ]: 47822681 : ABI_FREE(kpgx)
2984 [ - + ]: 47822681 : ABI_FREE(parity)
2985 [ - + ]: 47822681 : ABI_FREE(scalars)
2986 [ - + ]: 47822681 : ABI_FREE(teffv)
2987 : : !!$OMP END PARALLEL
2988 : :
2989 : 47822681 : end subroutine opernl4b
2990 : : !!***
2991 : :
2992 : : end module m_opernl
2993 : : !!***
|
