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 48329359 : subroutine opernl4a(choice,dgxdis,dgxds,d2gxdis,d2gxds2,dgxdt,&
1324 48329359 : & ffnl,gmet,gxa,ia3,idir,indlmn,ispinor,istwf_k,itypat,&
1325 48329359 : & jproj,kg_k,kpg_k,kpt,lmnmax,matblk,mincat,mlang1,mlang3,mlang4,&
1326 : & mlang5,mlang6,mproj,ndgxdt,nffnl,nincat,nkpg,nlang,nloalg,npw,&
1327 48329359 : & 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 48329359 : integer,allocatable :: parity(:)
1357 48329359 : 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 48329359 : 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 48329359 : jump=4
1370 :
1371 : !Get the actual maximum number of projectors
1372 341508396 : mmproj=maxval(indlmn(3,:,itypat))
1373 :
1374 : !Initialisation before blocking on the plane waves
1375 :
1376 : !Put projected scalars to zero
1377 10421732499 : gxa(:,:,:,1:mmproj)=0.0d0
1378 1893363449 : if (choice==2 .or. choice==4 .or. choice==5 .or. choice==23) dgxdt(:,:,:,:,1:mmproj)=0.0d0
1379 495190967 : if (choice==3 .or. choice==6 .or. choice==23) dgxds(:,:,:,1:mmproj)=0.0d0
1380 48329359 : 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 48329359 : ntens=1
1391 49553436 : if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5 .or. choice==23)ntens=4
1392 48329359 : if(nlang>=3 .or. (choice==3.or.choice==23))ntens=10
1393 48329359 : if(nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) )ntens=20
1394 48329359 : if(((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6)ntens=35
1395 48329359 : if(((choice==3.or.choice==23) .and. nlang==4) .or. (choice==6 .and. nlang>=2))ntens=56
1396 48329359 : if(choice==6 .and. nlang>=3)ntens=84
1397 48329359 : if(choice==6 .and. nlang==4)ntens=120
1398 :
1399 : !Set up second dimension of ffkg array, and allocate
1400 48329359 : nffkg=0 ; nffkge=0 ; nffkgd=0 ; nffkgk=0 ; nffkgs=0 ; nffkgs2=0
1401 166465836 : do ilang=1,nlang
1402 : ! Get the number of projectors for that angular momentum
1403 118136477 : 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 166465836 : if(nproj>0)then
1413 111662998 : ilang2=(ilang*(ilang+1))/2
1414 111662998 : nffkge=nffkge+nproj*ilang2
1415 111662998 : if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
1416 13948629 : 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 2883197 : ilang3=((ilang+2)*(ilang+3))/2
1419 2883197 : nffkgs=nffkgs+nproj*ilang3
1420 : end if
1421 111662998 : 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 48329359 : 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 193317436 : ABI_MALLOC(ffkg,(nffkg,mblkpw))
1451 144988077 : ABI_MALLOC(parity,(nffkg))
1452 144988077 : ABI_MALLOC(kpgx,(mblkpw,ntens))
1453 144988077 : ABI_MALLOC(scalars,(2,nffkg))
1454 48329359 : ABI_MALLOC(teffv,(2,mblkpw))
1455 :
1456 : !!$OMP DO
1457 140808733 : do ipw1=1,npw,mblkpw
1458 :
1459 92479374 : ipw2=min(npw,ipw1+mblkpw-1)
1460 92479374 : 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 92479374 : & npw,ntens,ntypat,parity)
1466 :
1467 287409519 : do ia=1,nincat
1468 :
1469 : ! Compute the shift eventually needed to get the phases in ph3d
1470 146600786 : iaph3d=ia
1471 146600786 : if(nloalg(2)>0)iaph3d=ia+ia3-1
1472 :
1473 2511622138 : do iffkg=1,nffkg
1474 2511622138 : 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 146600786 : 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 1792041487 : do ipw=1,nincpw
1510 1778981949 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1511 1778981949 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1512 1778981949 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1513 1792041487 : ig=ig+1
1514 : end do
1515 13059538 : 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 1580483160 : do ipw=1,nincpw
1558 1567777566 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1559 1567777566 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1560 1567777566 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1561 1567777566 : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1562 1567777566 : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1563 1567777566 : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1564 1580483160 : ig=ig+1
1565 : end do
1566 12705594 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1567 12705594 : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1568 12705594 : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1569 12705594 : 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 28134720 : 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 18537998016 : do ipw=1,nincpw
1751 18419531950 : ar=vect(1,ig)*ph3d(1,ig,iaph3d)-vect(2,ig)*ph3d(2,ig,iaph3d)
1752 18419531950 : ai=vect(2,ig)*ph3d(1,ig,iaph3d)+vect(1,ig)*ph3d(2,ig,iaph3d)
1753 18419531950 : teffv(1,ipw)=ar ; teffv(2,ipw)=ai
1754 18419531950 : scr1=scr1+ar*ffkg(1,ipw) ; sci1=sci1+ai*ffkg(1,ipw)
1755 18419531950 : scr2=scr2+ar*ffkg(2,ipw) ; sci2=sci2+ai*ffkg(2,ipw)
1756 18419531950 : scr3=scr3+ar*ffkg(3,ipw) ; sci3=sci3+ai*ffkg(3,ipw)
1757 18419531950 : scr4=scr4+ar*ffkg(4,ipw) ; sci4=sci4+ai*ffkg(4,ipw)
1758 18537998016 : ig=ig+1
1759 : end do
1760 118466066 : scalars(1,1)=scr1 ; scalars(2,1)=sci1
1761 118466066 : scalars(1,2)=scr2 ; scalars(2,2)=sci2
1762 118466066 : scalars(1,3)=scr3 ; scalars(2,3)=sci3
1763 118466066 : 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 597354970 : do start=1+jump,nffkg,jump
1888 478888904 : chunk=min(jump,nffkg-start+1)
1889 :
1890 118466066 : 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 3561863189 : do ipw=1,nincpw
1906 3539793021 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1907 3539793021 : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1908 3561863189 : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1909 : end do
1910 22070168 : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1911 22070168 : 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 65088399745 : do ipw=1,nincpw
1935 64650390007 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
1936 64650390007 : scr1=scr1+ar*ffkg(start ,ipw) ; sci1=sci1+ai*ffkg(start ,ipw)
1937 64650390007 : scr2=scr2+ar*ffkg(start+1,ipw) ; sci2=sci2+ai*ffkg(start+1,ipw)
1938 64650390007 : scr3=scr3+ar*ffkg(start+2,ipw) ; sci3=sci3+ai*ffkg(start+2,ipw)
1939 65088399745 : scr4=scr4+ar*ffkg(start+3,ipw) ; sci4=sci4+ai*ffkg(start+3,ipw)
1940 : end do
1941 438009738 : scalars(1,start )=scr1 ; scalars(2,start )=sci1
1942 438009738 : scalars(1,start+1)=scr2 ; scalars(2,start+1)=sci2
1943 438009738 : scalars(1,start+2)=scr3 ; scalars(2,start+2)=sci3
1944 438009738 : 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 478888904 : 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 146600786 : 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 277019944 : do iffkg=1,nffkg
2069 277019944 : scalars(parity(iffkg),iffkg)=0.0d0
2070 : end do
2071 : end if
2072 :
2073 146600786 : iffkg=0 ; iffkgs=nffkge+nffkgd ; iffkgk=nffkge*2
2074 146600786 : iffkgs2=nffkge+nffkgs
2075 600072299 : do ilang=1,nlang
2076 360992139 : nproj=jproj(ilang)
2077 507592925 : if(nproj>0)then
2078 : ! ilang2 is the number of independent tensor components
2079 : ! for symmetric tensor of rank ilang-1
2080 342798831 : ilang2=(ilang*(ilang+1))/2
2081 :
2082 : ! Loop over projectors
2083 986813988 : do iproj=1,nproj
2084 : ! Multiply by the k+G factors (tensors of various rank)
2085 2586630443 : do ii=1,ilang2
2086 : ! Get the starting address for the relevant tensor
2087 1942615286 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
2088 1942615286 : iffkg=iffkg+1
2089 : ! !$OMP CRITICAL (OPERNL4a_1)
2090 1942615286 : gxa(1,jj,ia,iproj)=gxa(1,jj,ia,iproj)+scalars(1,iffkg)
2091 1942615286 : 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 1942615286 : if ((choice==2.or.choice==23) .and. ndgxdt==3) then
2095 100848552 : do mu=1,3
2096 75636414 : jffkg=nffkge+(iffkg-1)*3+mu
2097 : ! Pay attention to the use of reals and imaginary parts here ...
2098 : ! !$OMP CRITICAL (OPERNL4a_2)
2099 75636414 : dgxdt(1,mu,jj,ia,iproj)=dgxdt(1,mu,jj,ia,iproj)-two_pi*scalars(2,jffkg)
2100 100848552 : 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 1942615286 : 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 2586630443 : 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 10554490 : ilang3=((ilang+2)*(ilang+3))/2
2138 10554490 : jjs=((ilang+1)*(ilang+2)*(ilang+3))/6
2139 : ! Compute strain gradient tensor components
2140 118974157 : do ii=1,ilang3
2141 : ! Note that iffkgs is also used by ddk and 2nd derivative parts
2142 108419667 : iffkgs=iffkgs+1
2143 108419667 : jj=ii+jjs
2144 : ! !$OMP CRITICAL (OPERNL4a_6)
2145 108419667 : dgxds(1,jj-4,ia,iproj)=dgxds(1,jj-4,ia,iproj)+scalars(1,iffkgs)
2146 118974157 : 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 644015157 : 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 986813988 : if (choice==5) then
2198 : ! Compute additional tensors related to ddk with ffnl(:,2,...)
2199 7759523 : ilangx=(ilang*(ilang+1))/2
2200 7759523 : jjs=((ilang-1)*ilang*(ilang+1))/6
2201 31704092 : do ii=1,ilangx
2202 : ! Note that iffkgs is also used by strain part
2203 23944569 : iffkgs=iffkgs+1
2204 23944569 : jj=ii+jjs
2205 : ! !$OMP CRITICAL (OPERNL4a_10)
2206 23944569 : dgxdt(1,1,jj,ia,iproj)=dgxdt(1,1,jj,ia,iproj)+scalars(1,iffkgs)
2207 31704092 : 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 7759523 : if(ilang>=2)then
2212 4953587 : ilangx=((ilang-1)*ilang)/2
2213 4953587 : jjs=((ilang-2)*(ilang-1)*ilang)/6
2214 14093562 : do ii=1,ilangx
2215 9139975 : iffkgk=iffkgk+1
2216 9139975 : jj=ii+jjs
2217 : ! !$OMP CRITICAL (OPERNL4a_11)
2218 9139975 : dgxdt(1,2,jj,ia,iproj)=dgxdt(1,2,jj,ia,iproj)+scalars(1,iffkgk)
2219 14093562 : 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 48329359 : ABI_FREE(ffkg)
2242 48329359 : ABI_FREE(kpgx)
2243 48329359 : ABI_FREE(parity)
2244 48329359 : ABI_FREE(scalars)
2245 48329359 : ABI_FREE(teffv)
2246 : !!$OMP END PARALLEL
2247 :
2248 : !DEBUG
2249 : !write(std_out,*)' opernl4a : exit'
2250 : !ENDDEBUG
2251 :
2252 48329359 : 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 47642811 : subroutine opernl4b(choice,dgxds,dgxdt,ffnl,gmet,gxa,&
2326 47642811 : & ia3,idir,indlmn,ispinor,itypat,jproj,kg_k,kpg_k,kpt,&
2327 : & lmnmax,matblk,mincat,mlang3,mlang4,mproj,ndgxdt,nffnl,nincat,&
2328 47642811 : & 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 47642811 : integer,allocatable :: parity(:)
2356 47642811 : real(dp),allocatable :: ffkg(:,:),kpgx(:,:),scalars(:,:),teffv(:,:)
2357 :
2358 : ! *************************************************************************
2359 :
2360 : !mblkpw sets the size of blocks of planewaves
2361 47642811 : 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 47642811 : 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 47642811 : ntens=1
2375 47642811 : if(nlang>=2 .or. choice==2 .or. choice==4 .or. choice==5)ntens=4
2376 47642811 : if(nlang>=3 .or. choice==3)ntens=10
2377 47642811 : if(nlang>=4 .or. (choice==3 .and. nlang>=2) )ntens=20
2378 47642811 : if(choice==3 .and. nlang>=3)ntens=35
2379 47642811 : if(choice==3 .and. nlang==4)ntens=56
2380 :
2381 : !Set up second dimension of ffkg array, and allocate
2382 47642811 : nffkg=0; nffkge=0; nffkgd=0; nffkgk=0; nffkgs=0
2383 163979918 : do ilang=1,nlang
2384 : ! Get the number of projectors for that angular momentum
2385 116337107 : nproj=jproj(ilang)
2386 : ! If there is a non-local part, accumulate the number of vectors needed
2387 163979918 : if(nproj>0)then
2388 109923847 : ilang2=(ilang*(ilang+1))/2
2389 109923847 : nffkge=nffkge+nproj*ilang2
2390 109923847 : if(choice==5)nffkgk=nffkgk+nproj*(2*ilang2-ilang)
2391 109923847 : if(choice==2 .or. choice==4)nffkgd=nffkgd+ndgxdt*nproj*ilang2
2392 109923847 : 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 47642811 : 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 190571244 : ABI_MALLOC(ffkg,(nffkg,mblkpw))
2407 142928433 : ABI_MALLOC(parity,(nffkg))
2408 142928433 : ABI_MALLOC(kpgx,(mblkpw,ntens))
2409 142928433 : ABI_MALLOC(scalars,(2,nffkg))
2410 47642811 : ABI_MALLOC(teffv,(2,mblkpw))
2411 :
2412 : !Loop on subsets of plane waves (blocking)
2413 : !!$OMP DO
2414 138729023 : do ipw1=1,npw,mblkpw
2415 :
2416 91086212 : ipw2=min(npw,ipw1+mblkpw-1)
2417 91086212 : 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 91086212 : & npw,ntens,ntypat,parity)
2423 :
2424 138729023 : 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 233428345 : do ia=1,nincat
2432 :
2433 : ! Compute the shift eventually needed to get the phases in ph3d
2434 142342133 : iaph3d=ia
2435 142342133 : if(nloalg(2)>0)iaph3d=ia+ia3-1
2436 :
2437 : ! Transfer gxa (and eventually dgxdt) in scalars with different indexing
2438 142342133 : iffkg=0
2439 142342133 : iffkgk=nffkge*2
2440 142342133 : iffkgs=nffkge
2441 492226992 : do ilang=1,nlang
2442 349884859 : nproj=jproj(ilang)
2443 492226992 : if (nproj>0) then
2444 331948266 : ilang2=(ilang*(ilang+1))/2
2445 331948266 : ilang3=((ilang+2)*(ilang+3))/2
2446 955133184 : do iproj=1,nproj
2447 2499476602 : do ii=1,ilang2
2448 1876291684 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
2449 1876291684 : iffkg=iffkg+1
2450 2499476602 : if(choice==1 .or. choice==3)then
2451 1707903956 : scalars(1,iffkg)=gxa(1,jj,ia,iproj)
2452 1707903956 : scalars(2,iffkg)=gxa(2,jj,ia,iproj)
2453 168387728 : 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 23944569 : else if (choice==5) then
2464 23944569 : jffkg=nffkge+iffkg
2465 : ! The gxa and dgxdt arrays are switched, in order
2466 : ! to give the correct combination when multiplying ffkg,
2467 23944569 : scalars(1,jffkg)= gxa(1,jj,ia,iproj)
2468 23944569 : scalars(2,jffkg)= gxa(2,jj,ia,iproj)
2469 23944569 : scalars(1,iffkg)= dgxdt(1,1,jj,ia,iproj)
2470 23944569 : scalars(2,iffkg)= dgxdt(2,1,jj,ia,iproj)
2471 : end if
2472 : end do
2473 623184918 : 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 955133184 : if(ilang>=2 .and. choice==5)then
2482 14093562 : do ii=1,((ilang-1)*ilang)/2
2483 9139975 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
2484 9139975 : iffkgk=iffkgk+1
2485 9139975 : scalars(1,iffkgk)= dgxdt(1,2,jj,ia,iproj)
2486 14093562 : 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 233428345 : if(nffkg<=jump)then
2509 :
2510 13054890 : select case(nffkg)
2511 :
2512 : case(1)
2513 :
2514 13054890 : scr1=scalars(1,1) ; sci1=scalars(2,1)
2515 13054890 : ig=ipw1
2516 1791422219 : do ipw=1,nincpw
2517 1778367329 : ar=ffkg(1,ipw)*scr1 ; ai=ffkg(1,ipw)*sci1
2518 1778367329 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2519 1778367329 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2520 1791422219 : 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 12627080 : scr1=scalars(1,1) ; sci1=scalars(2,1)
2554 12627080 : scr2=scalars(1,2) ; sci2=scalars(2,2)
2555 12627080 : scr3=scalars(1,3) ; sci3=scalars(2,3)
2556 12627080 : scr4=scalars(1,4) ; sci4=scalars(2,4)
2557 12627080 : ig=ipw1
2558 1566526631 : do ipw=1,nincpw
2559 1553899551 : ar= ffkg(1,ipw)*scr1 ; ai= ffkg(1,ipw)*sci1
2560 1553899551 : ar=ar+ffkg(2,ipw)*scr2 ; ai=ai+ffkg(2,ipw)*sci2
2561 1553899551 : ar=ar+ffkg(3,ipw)*scr3 ; ai=ai+ffkg(3,ipw)*sci3
2562 1553899551 : ar=ar+ffkg(4,ipw)*scr4 ; ai=ai+ffkg(4,ipw)*sci4
2563 1553899551 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2564 1553899551 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2565 1566526631 : 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 28051558 : 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 17928306203 : do ipw=1,nincpw
2662 17928306203 : teffv(1,ipw)=0.0d0 ; teffv(2,ipw)=0.0d0
2663 : end do
2664 :
2665 : ! 2) Here treart the loop for big nffkg
2666 532563899 : do start=1,nffkg-jump,jump
2667 418273324 : chunk=min(jump,nffkg-jump-start+1)
2668 :
2669 114290575 : 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 21257686 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2683 21257686 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2684 3440293509 : do ipw=1,nincpw
2685 3419035823 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2686 3419035823 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2687 3419035823 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2688 3440293509 : 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 378411471 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2707 378411471 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2708 378411471 : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2709 378411471 : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2710 56247294316 : do ipw=1,nincpw
2711 55868882845 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2712 55868882845 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2713 55868882845 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2714 55868882845 : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2715 55868882845 : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2716 56247294316 : 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 418273324 : 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 114290575 : 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 114290575 : scr1=scalars(1,start ) ; sci1=scalars(2,start )
2854 114290575 : scr2=scalars(1,start+1) ; sci2=scalars(2,start+1)
2855 114290575 : scr3=scalars(1,start+2) ; sci3=scalars(2,start+2)
2856 114290575 : scr4=scalars(1,start+3) ; sci4=scalars(2,start+3)
2857 114290575 : ig=ipw1
2858 17928306203 : do ipw=1,nincpw
2859 17814015628 : ar=teffv(1,ipw) ; ai=teffv(2,ipw)
2860 17814015628 : ar=ar+ffkg(start ,ipw)*scr1 ; ai=ai+ffkg(start ,ipw)*sci1
2861 17814015628 : ar=ar+ffkg(start+1,ipw)*scr2 ; ai=ai+ffkg(start+1,ipw)*sci2
2862 17814015628 : ar=ar+ffkg(start+2,ipw)*scr3 ; ai=ai+ffkg(start+2,ipw)*sci3
2863 17814015628 : ar=ar+ffkg(start+3,ipw)*scr4 ; ai=ai+ffkg(start+3,ipw)*sci4
2864 17814015628 : vect(1,ig)=vect(1,ig)+ar*ph3d(1,ig,iaph3d)+ai*ph3d(2,ig,iaph3d)
2865 17814015628 : vect(2,ig)=vect(2,ig)+ai*ph3d(1,ig,iaph3d)-ar*ph3d(2,ig,iaph3d)
2866 17928306203 : 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 47642811 : ABI_FREE(ffkg)
2983 47642811 : ABI_FREE(kpgx)
2984 47642811 : ABI_FREE(parity)
2985 47642811 : ABI_FREE(scalars)
2986 47642811 : ABI_FREE(teffv)
2987 : !!$OMP END PARALLEL
2988 :
2989 47642811 : end subroutine opernl4b
2990 : !!***
2991 :
2992 : end module m_opernl
2993 : !!***
|
