Line data Source code
1 : !!****m* ABINIT/m_opernlc_ylm_allwf
2 : !! NAME
3 : !! m_opernlc_ylm_allwf
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2008-2026 ABINIT group (MT)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_opernlc_ylm_allwf
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 : use m_xmpi
27 : use m_gputk
28 : use m_abi_linalg
29 : use, intrinsic :: iso_c_binding
30 :
31 : use defs_abitypes, only : MPI_type
32 : use m_opernlc_ylm, only : ls_ylm
33 :
34 : implicit none
35 :
36 : private
37 : !!***
38 :
39 : public :: opernlc_ylm_allwf
40 : !!***
41 :
42 : ! Work buffers to be used when iphase==2
43 : real(dp), allocatable, target :: d2gxdtfac_2ndphase(:,:,:,:,:)
44 : real(dp), allocatable, target :: dgxdtfac_2ndphase(:,:,:,:,:)
45 : real(dp), allocatable, target :: gxfac_2ndphase(:,:,:,:)
46 : ! Work buffer for NC+SO L.S matrix (pointer pattern for compiler robustness)
47 : real(dp), allocatable, target :: ls_ylm_so_data(:,:,:)
48 :
49 : !----------------------------------------------------------------------
50 :
51 : contains
52 : !!***
53 :
54 : !----------------------------------------------------------------------
55 :
56 : !!****f* m_opernlc_ylm_allwf/alloc_work_arrays
57 : !! NAME
58 : !! alloc_work_arrays
59 : !!
60 : !! FUNCTION
61 : !! Allocation of work arrays
62 : !!
63 : !! INPUTS
64 : !!
65 : !! SOURCE
66 0 : subroutine alloc_work_arrays(optder,cplex_fac,ndgxdtfac,nd2gxdtfac,nprojs,nspinor,ndat,gpu_option)
67 :
68 : integer,intent(in) :: optder,cplex_fac,ndgxdtfac,nd2gxdtfac,nprojs,nspinor,ndat,gpu_option
69 :
70 : ! *************************************************************************
71 :
72 0 : ABI_MALLOC(gxfac_2ndphase,(cplex_fac,nprojs,nspinor,ndat))
73 : #ifdef HAVE_OPENMP_OFFLOAD
74 : !$OMP TARGET ENTER DATA MAP(alloc:gxfac_2ndphase) IF(gpu_option==ABI_GPU_OPENMP)
75 : #endif
76 0 : if(gpu_option==ABI_GPU_OPENMP) then
77 0 : call gpu_set_to_zero(gxfac_2ndphase, int(cplex_fac,c_size_t)*nprojs*nspinor*ndat)
78 : else
79 0 : gxfac_2ndphase(:,:,:,:) = zero
80 : end if
81 0 : if (optder>=1) then
82 0 : ABI_MALLOC(dgxdtfac_2ndphase,(cplex_fac,ndgxdtfac,nprojs,nspinor,ndat))
83 : #ifdef HAVE_OPENMP_OFFLOAD
84 : !$OMP TARGET ENTER DATA MAP(alloc:dgxdtfac_2ndphase) IF(gpu_option==ABI_GPU_OPENMP)
85 : #endif
86 0 : if(gpu_option==ABI_GPU_OPENMP) then
87 0 : call gpu_set_to_zero(dgxdtfac_2ndphase, int(cplex_fac,c_size_t)*ndgxdtfac*nprojs*nspinor*ndat)
88 : else
89 0 : dgxdtfac_2ndphase(:,:,:,:,:) = zero
90 : end if
91 : end if
92 0 : if (optder>=2) then
93 0 : ABI_MALLOC(d2gxdtfac_2ndphase,(cplex_fac,nd2gxdtfac,nprojs,nspinor,ndat))
94 : #ifdef HAVE_OPENMP_OFFLOAD
95 : !$OMP TARGET ENTER DATA MAP(alloc:dgxdtfac_2ndphase) IF(gpu_option==ABI_GPU_OPENMP)
96 : #endif
97 0 : if(gpu_option==ABI_GPU_OPENMP) then
98 0 : call gpu_set_to_zero(d2gxdtfac_2ndphase, int(cplex_fac,c_size_t)*nd2gxdtfac*nprojs*nspinor*ndat)
99 : else
100 0 : d2gxdtfac_2ndphase(:,:,:,:,:) = zero
101 : end if
102 : end if
103 :
104 0 : end subroutine alloc_work_arrays
105 : !!***
106 :
107 : !----------------------------------------------------------------------
108 :
109 : !!****f* m_opernlc_ylm_allwf/destroy_work_arrays
110 : !! NAME
111 : !! destroy_work_arrays
112 : !!
113 : !! FUNCTION
114 : !! Destruction of work arrays
115 : !!
116 : !! INPUTS
117 : !!
118 : !! SOURCE
119 0 : subroutine destroy_work_arrays(gpu_option)
120 :
121 : integer,intent(in) :: gpu_option
122 :
123 : ! *************************************************************************
124 :
125 : ABI_UNUSED(gpu_option) !Silent abirules
126 :
127 0 : if(allocated(gxfac_2ndphase)) then
128 : #ifdef HAVE_OPENMP_OFFLOAD
129 : !$OMP TARGET EXIT DATA MAP(delete:gxfac_2ndphase) IF(gpu_option==ABI_GPU_OPENMP)
130 : #endif
131 0 : ABI_FREE(gxfac_2ndphase)
132 : end if
133 0 : if(allocated(dgxdtfac_2ndphase)) then
134 : #ifdef HAVE_OPENMP_OFFLOAD
135 : !$OMP TARGET EXIT DATA MAP(delete:dgxdtfac_2ndphase) IF(gpu_option==ABI_GPU_OPENMP)
136 : #endif
137 0 : ABI_FREE(dgxdtfac_2ndphase)
138 : end if
139 0 : if(allocated(d2gxdtfac_2ndphase)) then
140 : #ifdef HAVE_OPENMP_OFFLOAD
141 : !$OMP TARGET EXIT DATA MAP(delete:d2gxdtfac_2ndphase) IF(gpu_option==ABI_GPU_OPENMP)
142 : #endif
143 0 : ABI_FREE(d2gxdtfac_2ndphase)
144 : end if
145 :
146 0 : end subroutine destroy_work_arrays
147 : !!***
148 :
149 : !----------------------------------------------------------------------
150 :
151 : !!****f* ABINIT/opernlc_ylm_allwf
152 : !! NAME
153 : !! opernlc_ylm_allwf
154 : !!
155 : !! FUNCTION
156 : !! * Operate with the non-local part of the hamiltonian,
157 : !! in order to reduce projected scalars
158 : !! * Operate with the non-local projectors and the overlap matrix,
159 : !! in order to reduce projected scalars
160 : !!
161 : !! INPUTS
162 : !! atindx1(natom)=index table for atoms (gives the absolute index of
163 : !! an atom from its rank in a block of atoms)
164 : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
165 : !! 2 if <p_lmn|c> scalars are complex
166 : !! cplex_dgxdt(ndgxdt) = used only when cplex = 1
167 : !! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
168 : !! cplex_enl=1 if enl factors are real, 2 if they are complex
169 : !! cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
170 : !! dgxdt(cplex,ndgxdt,nlmn,nincat)=grads of projected scalars (only if optder>0)
171 : !! dimenl1,dimenl2=dimensions of enl (see enl)
172 : !! dimekbq=1 if enl factors do not contain a exp(-iqR) phase, 2 is they do
173 : !! enl(cplex_enl*dimenl1,dimenl2,nspinortot**2,dimekbq)=
174 : !! ->Norm conserving : ==== when paw_opt=0 ====
175 : !! (Real) Kleinman-Bylander energies (hartree)
176 : !! dimenl1=lmnmax - dimenl2=ntypat
177 : !! dimekbq is 2 if Enl contains a exp(-iqR) phase, 1 otherwise
178 : !! ->PAW : ==== when paw_opt=1, 2 or 4 ====
179 : !! (Real or complex, hermitian) Dij coefs to connect projectors
180 : !! dimenl1=cplex_enl*lmnmax*(lmnmax+1)/2 - dimenl2=natom
181 : !! These are complex numbers if cplex_enl=2
182 : !! enl(:,:,1) contains Dij^up-up
183 : !! enl(:,:,2) contains Dij^dn-dn
184 : !! enl(:,:,3) contains Dij^up-dn (only if nspinor=2)
185 : !! enl(:,:,4) contains Dij^dn-up (only if nspinor=2)
186 : !! dimekbq is 2 if Dij contains a exp(-iqR) phase, 1 otherwise
187 : !! gx(cplex,nlmn,nincat*abs(enl_opt))= projected scalars
188 : !! iatm=absolute rank of first atom of the current block of atoms
189 : !! indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
190 : !! itypat=type of atoms
191 : !! lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
192 : !! mpi_enreg=information about MPI parallelization
193 : !! natom=number of atoms in cell
194 : !! ndgxdt=second dimension of dgxdt
195 : !! ndgxdtfac=second dimension of dgxdtfac
196 : !! nincat=number of atoms in the subset here treated
197 : !! nlmn=number of (l,m,n) numbers for current type of atom
198 : !! nspinor= number of spinorial components of the wavefunctions (on current proc)
199 : !! nspinortot=total number of spinorial components of the wavefunctions
200 : !! optder=0=only gxfac is computed, 1=both gxfac and dgxdtfac are computed
201 : !! 2=gxfac, dgxdtfac and d2gxdtfac are computed
202 : !! paw_opt= define the nonlocal operator concerned with:
203 : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
204 : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
205 : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
206 : !! paw_opt=3 : PAW overlap matrix (Sij)
207 : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
208 : !! sij(nlm*(nlmn+1)/2)=overlap matrix components (only if paw_opt=2, 3 or 4)
209 : !!
210 : !! OUTPUT
211 : !! if (paw_opt=0, 1, 2 or 4)
212 : !! gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
213 : !! if (paw_opt=3 or 4)
214 : !! gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
215 : !! if (optder==1.and.paw_opt=0, 1, 2 or 4)
216 : !! dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Vnl (NL operator)
217 : !! if (optder==1.and.paw_opt=3 or 4)
218 : !! dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Sij (overlap)
219 : !!
220 : !! NOTES
221 : !! This routine operates for one type of atom, and within this given type of atom,
222 : !! for a subset of at most nincat atoms.
223 : !!
224 : !! About the non-local factors symmetry:
225 : !! - The lower triangular part of the Dij matrix can be deduced from the upper one
226 : !! with the following relation: D^s2s1_ji = (D^s1s2_ij)^*
227 : !! where s1,s2 are spinor components
228 : !! - The Dij factors can contain a exp(-iqR) phase
229 : !! This phase does not have to be included in the symmetry rule
230 : !! For that reason, we first apply the real part (cos(qR).D^s1s2_ij)
231 : !! then, we apply the imaginary part (-sin(qR).D^s1s2_ij)
232 : !!
233 : !! SOURCE
234 :
235 36992 : subroutine opernlc_ylm_allwf(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,&
236 36992 : & dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,dimekbq,enl,&
237 36992 : & gx,gxfac,gxfac_sij,iatm,indlmn,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,&
238 18496 : & nd2gxdt,nd2gxdtfac,nincat,nlmn,nspinor,nspinortot,optder,paw_opt,sij,ndat,ibeg,iend,nprojs,ndat_enl,gpu_option)
239 :
240 : !Arguments ------------------------------------
241 : !scalars
242 : integer,intent(in) :: cplex,cplex_enl,cplex_fac,dimenl1,dimenl2,dimekbq,iatm,itypat
243 : integer,intent(in) :: natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,nincat,nspinor,nspinortot,optder,paw_opt,gpu_option
244 : integer,intent(inout) :: nlmn
245 : integer,intent(in) :: ndat,ibeg,iend,nprojs,ndat_enl
246 : real(dp) :: lambda(ndat)
247 : type(MPI_type) , intent(in) :: mpi_enreg
248 : !arrays
249 : integer,intent(in) :: atindx1(natom),indlmn(6,nlmn),cplex_dgxdt(ndgxdt),cplex_d2gxdt(nd2gxdt)
250 : real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nprojs,nspinor,ndat)
251 : real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor,ndat)
252 : real(dp),intent(in),target :: enl(dimenl1,dimenl2,nspinortot**2,ndat_enl,dimekbq)
253 : real(dp),intent(inout) :: gx(cplex,nprojs,nspinor,ndat)
254 : real(dp),intent(in) :: sij(:)
255 : real(dp),intent(out),target :: dgxdtfac(cplex_fac,ndgxdtfac,nprojs,nspinor,ndat)
256 : real(dp),intent(out) :: dgxdtfac_sij(cplex,ndgxdtfac,nprojs,nspinor,ndat*(paw_opt/3))
257 : real(dp),intent(out),target :: d2gxdtfac(cplex_fac,nd2gxdtfac,nprojs,nspinor,ndat)
258 : real(dp),intent(out) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nprojs,nspinor,ndat*(paw_opt/3))
259 : real(dp),intent(out),target :: gxfac(cplex_fac,nprojs,nspinor,ndat)
260 : real(dp),intent(out) :: gxfac_sij(cplex,nprojs,nspinor,ndat)
261 :
262 : !Local variables-------------------------------
263 : !Arrays
264 : !scalars
265 : integer :: cplex_,ia,ijlmn,ilm,ilmn,i0lmn,iln,index_enl,iphase,ispinor,ispinor_index,idat
266 : integer :: jlm,j0lmn,jjlmn,jlmn,jspinor,mu,shift,ii
267 : integer :: ll_so,klm_so,lmax_so,nlmso,sign_so
268 : real(dp) :: ekb_so_1,ekb_so_2,ls_uu_im,ls_ud_re,ls_ud_im
269 : !arrays
270 36992 : real(dp) :: enl_(2),gxfi(2),gxi(cplex),gxj(cplex)
271 18496 : real(dp), ABI_CONTIGUOUS pointer :: d2gxdtfac_(:,:,:,:,:),dgxdtfac_(:,:,:,:,:),gxfac_(:,:,:,:)
272 18496 : real(dp), ABI_CONTIGUOUS pointer :: ls_ylm_so_(:,:,:)
273 18496 : real(dp), ABI_CONTIGUOUS pointer :: enl_ptr(:,:,:),enl_ptr2(:,:,:,:)
274 :
275 : ! *************************************************************************
276 :
277 0 : if (gpu_option/=ABI_GPU_DISABLED.and.mpi_enreg%paral_spinor==1) then
278 0 : ABI_ERROR('parallelization over spinors (npspinor=2) not allowed with GPU!')
279 : end if
280 :
281 : ABI_UNUSED(iend)
282 : ABI_UNUSED(d2gxdt)
283 : ABI_UNUSED(cplex_d2gxdt)
284 : ABI_UNUSED(d2gxdtfac_sij)
285 : DBG_ENTER("COLL")
286 :
287 : !Parallelization over spinors treatment
288 18496 : shift=0;if (mpi_enreg%paral_spinor==1) shift=mpi_enreg%me_spinor
289 : !When Enl factors contain a exp(-iqR) phase:
290 : ! - We loop over the real and imaginary parts
291 : ! - We need an additional memory space
292 36992 : do iphase=1,dimekbq
293 18496 : if (paw_opt==3) cycle
294 18496 : if (iphase==1) then
295 18496 : gxfac_ => gxfac ; dgxdtfac_ => dgxdtfac ; d2gxdtfac_ => d2gxdtfac
296 : else
297 0 : ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac==1 when dimekbq=2!")
298 0 : call alloc_work_arrays(optder,cplex_fac,ndgxdtfac,nd2gxdtfac,nprojs,nspinor,ndat,gpu_option)
299 0 : gxfac_ => gxfac_2ndphase
300 0 : if(optder>=1) dgxdtfac_ => dgxdtfac_2ndphase
301 0 : if(optder>=2) d2gxdtfac_ => d2gxdtfac_2ndphase
302 : end if
303 18496 : enl_ptr => enl(:,:,:,1,iphase)
304 18496 : enl_ptr2 => enl(:,:,:,:,iphase)
305 :
306 : !NC+SO: precompute L.S matrix once
307 18496 : lmax_so = 0
308 18496 : if (paw_opt==0.and.nspinortot==2.and.nspinor==nspinortot) then
309 0 : if (any(indlmn(6,1:nlmn)==2)) then
310 0 : do ilmn=1,nlmn
311 0 : if (indlmn(6,ilmn)==2) lmax_so = max(lmax_so, indlmn(1,ilmn))
312 : end do
313 0 : if (lmax_so > 0) then
314 0 : nlmso = (lmax_so+1)**2*((lmax_so+1)**2+1)/2
315 0 : ABI_MALLOC(ls_ylm_so_data,(2,nlmso,2))
316 0 : call ls_ylm(ls_ylm_so_data, lmax_so)
317 : #ifdef HAVE_OPENMP_OFFLOAD
318 : !$OMP TARGET ENTER DATA MAP(to:ls_ylm_so_data) IF(gpu_option==ABI_GPU_OPENMP)
319 : #endif
320 0 : ls_ylm_so_ => ls_ylm_so_data
321 : end if
322 : end if
323 : end if
324 18496 : if (paw_opt==0.and.mpi_enreg%paral_spinor==1.and.lmax_so>0) then
325 0 : ABI_ERROR('parallelization over spinors, spin-orbit and norm-conserving psps not implemented!')
326 : end if
327 :
328 :
329 : !Accumulate gxfac related to non-local operator (Norm-conserving)
330 : !-------------------------------------------------------------------
331 18496 : if (paw_opt==0) then
332 :
333 : !Enl is E(Kleinman-Bylander)
334 3200 : ABI_CHECK(cplex_enl/=2,"BUG: invalid cplex_enl=2!")
335 3200 : ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!")
336 :
337 3200 : if (lmax_so == 0) then
338 :
339 : ! NC+SR ---
340 : #ifdef HAVE_OPENMP_OFFLOAD
341 : !$OMP TARGET TEAMS DISTRIBUTE &
342 : !$OMP& MAP(to:gxfac_,gx,enl_ptr2,indlmn) &
343 : !$OMP& PRIVATE(idat,ispinor) &
344 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
345 : #endif
346 39552 : do idat=1,ndat
347 75904 : do ispinor=1,nspinor
348 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(ia,ilmn,iln,ii)
349 236288 : do ia=1,nincat
350 2599168 : do ilmn=1,nlmn
351 6761472 : do ii=1,cplex
352 6597888 : if (indlmn(6,ilmn)==2) then
353 0 : gxfac_(ii,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat) = zero
354 : else
355 4198656 : iln = indlmn(5,ilmn)
356 : gxfac_(ii,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
357 4198656 : & enl_ptr2(iln,itypat,ispinor+shift,min(ndat_enl,idat))*gx(ii,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
358 : end if
359 : end do
360 : end do
361 : end do
362 : end do
363 : end do
364 :
365 : else
366 :
367 : ! NC+SR+SO: real-Ylm L.S coupling ---
368 : #ifdef HAVE_OPENMP_OFFLOAD
369 : !$OMP TARGET TEAMS DISTRIBUTE &
370 : !$OMP& MAP(to:gxfac_,gx,enl_ptr2,indlmn,ls_ylm_so_) &
371 : !$OMP& PRIVATE(idat) &
372 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
373 : #endif
374 0 : do idat=1,ndat
375 : !$OMP PARALLEL DO COLLAPSE(2) &
376 : !$OMP& PRIVATE(ia,ilmn,iln,ekb_so_1,ekb_so_2,ll_so,ilm,jlmn,jlm,klm_so,sign_so,ls_uu_im,ls_ud_re,ls_ud_im)
377 0 : do ia=1,nincat
378 0 : do ilmn=1,nlmn
379 0 : iln = indlmn(5,ilmn)
380 :
381 0 : if (indlmn(6,ilmn)/=2) then
382 0 : ekb_so_1 = enl_ptr2(iln,itypat,1,min(ndat_enl,idat))
383 0 : ekb_so_2 = enl_ptr2(iln,itypat,2,min(ndat_enl,idat))
384 0 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat)=ekb_so_1*gx(1,ilmn+(ia-1)*nlmn+ibeg,1,idat)
385 0 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat)=ekb_so_1*gx(2,ilmn+(ia-1)*nlmn+ibeg,1,idat)
386 0 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,2,idat)=ekb_so_2*gx(1,ilmn+(ia-1)*nlmn+ibeg,2,idat)
387 0 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,2,idat)=ekb_so_2*gx(2,ilmn+(ia-1)*nlmn+ibeg,2,idat)
388 :
389 : else
390 0 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat) = zero
391 0 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat) = zero
392 0 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,2,idat) = zero
393 0 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,2,idat) = zero
394 :
395 0 : ekb_so_1 = enl_ptr2(iln,itypat,1,min(ndat_enl,idat))
396 0 : if (abs(ekb_so_1)<tol16) cycle
397 0 : ll_so = indlmn(1,ilmn)
398 0 : ilm = indlmn(4,ilmn)
399 0 : do jlmn=1,nlmn
400 0 : if (indlmn(6,jlmn)/=2) cycle
401 0 : if (indlmn(1,jlmn)/=ll_so) cycle
402 0 : if (indlmn(3,jlmn)/=indlmn(3,ilmn)) cycle
403 0 : jlm = indlmn(4,jlmn)
404 0 : if (ilm<=jlm) then
405 0 : klm_so = jlm*(jlm-1)/2 + ilm
406 0 : sign_so = 1
407 : else
408 0 : klm_so = ilm*(ilm-1)/2 + jlm
409 0 : sign_so = -1
410 : end if
411 0 : ls_uu_im = sign_so * ls_ylm_so_(2,klm_so,1)
412 0 : ls_ud_re = sign_so * ls_ylm_so_(1,klm_so,2)
413 0 : ls_ud_im = sign_so * ls_ylm_so_(2,klm_so,2)
414 :
415 : ! up-up
416 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
417 0 : & - ekb_so_1*ls_uu_im*gx(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)
418 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
419 0 : & + ekb_so_1*ls_uu_im*gx(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)
420 : ! up-dn
421 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
422 0 : & + ekb_so_1*(ls_ud_re*gx(1,jlmn+(ia-1)*nlmn+ibeg,2,idat) - ls_ud_im*gx(2,jlmn+(ia-1)*nlmn+ibeg,2,idat))
423 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
424 0 : & + ekb_so_1*(ls_ud_re*gx(2,jlmn+(ia-1)*nlmn+ibeg,2,idat) + ls_ud_im*gx(1,jlmn+(ia-1)*nlmn+ibeg,2,idat))
425 : ! dn-up
426 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,2,idat)=gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
427 0 : & + ekb_so_1*(-ls_ud_re*gx(1,jlmn+(ia-1)*nlmn+ibeg,1,idat) - ls_ud_im*gx(2,jlmn+(ia-1)*nlmn+ibeg,1,idat))
428 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,2,idat)=gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
429 0 : & + ekb_so_1*(-ls_ud_re*gx(2,jlmn+(ia-1)*nlmn+ibeg,1,idat) + ls_ud_im*gx(1,jlmn+(ia-1)*nlmn+ibeg,1,idat))
430 : ! dn-dn
431 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,2,idat)=gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
432 0 : & + ekb_so_1*ls_uu_im*gx(2,jlmn+(ia-1)*nlmn+ibeg,2,idat)
433 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,2,idat)=gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
434 0 : & - ekb_so_1*ls_uu_im*gx(1,jlmn+(ia-1)*nlmn+ibeg,2,idat)
435 : end do ! jlmn
436 : end if ! indlmn(:,6)==2
437 : end do ! ilmn
438 : end do ! ia
439 : !$OMP END PARALLEL DO
440 :
441 : end do ! idat
442 : end if ! NC+SO
443 : end if ! NC
444 :
445 : !Accumulate gxfac related to nonlocal operator (PAW)
446 : !-------------------------------------------------------------------
447 18496 : if (paw_opt==1.or.paw_opt==2.or.paw_opt==4) then
448 : !Enl is psp strength Dij or (Dij-lambda.Sij)
449 :
450 : ! === Diagonal term(s) (up-up, down-down)
451 :
452 : ! 1-Enl is real
453 15296 : if (cplex_enl==1) then
454 15296 : if (paw_opt==2) then
455 : #ifdef HAVE_OPENMP_OFFLOAD
456 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
457 : !$OMP& MAP(to:gxfac_,enl_ptr2,atindx1,gx,sij,lambda) &
458 : !$OMP& PRIVATE(idat,ispinor) &
459 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
460 : #endif
461 1576 : do idat=1,ndat
462 2568 : do ispinor=1,nspinor
463 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,ispinor_index,index_enl,j0lmn,i0lmn,ii)
464 3520 : do ia=1,nincat
465 19296 : do jlmn=1,nlmn
466 16768 : ispinor_index=ispinor+shift
467 16768 : index_enl=atindx1(iatm+ia)
468 16768 : j0lmn=jlmn*(jlmn-1)/2
469 46256 : do ii=1,cplex
470 : gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
471 : & gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat) + &
472 : & (enl_ptr2(j0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat))-lambda(idat) * sij(j0lmn+jlmn)) * &
473 46256 : & gx(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
474 : end do
475 115776 : do ilmn=1,jlmn-1
476 284504 : do ii=1,cplex
477 : gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
478 : & gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat) + &
479 : & (enl_ptr2(j0lmn+ilmn,index_enl,ispinor_index,min(ndat_enl,idat))-lambda(idat) * sij(j0lmn+ilmn)) * &
480 267736 : & gx(ii,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
481 : end do
482 : end do
483 18304 : if(jlmn<nlmn) then
484 114240 : do ilmn=jlmn+1,nlmn
485 99008 : i0lmn=(ilmn*(ilmn-1)/2)
486 282968 : do ii=1,cplex
487 : gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
488 : & gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat) + &
489 : & (enl_ptr2(i0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat))-lambda(idat) * sij(i0lmn+jlmn)) *&
490 267736 : & gx(ii,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
491 : end do
492 : end do
493 : end if
494 : end do
495 : end do
496 : end do
497 : end do
498 :
499 : else
500 : #ifdef HAVE_OPENMP_OFFLOAD
501 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
502 : !$OMP& MAP(to:enl_ptr2,atindx1,gx,gxfac_) &
503 : !$OMP& PRIVATE(idat,ispinor) &
504 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
505 : #endif
506 47206 : do idat=1,ndat
507 79700 : do ispinor=1,nspinor
508 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(jlmn,j0lmn,ii,ia,ispinor_index,index_enl)
509 107874 : do ia=1,nincat
510 639488 : do jlmn=1,nlmn
511 564108 : ispinor_index=ispinor+shift
512 564108 : index_enl=atindx1(iatm+ia)
513 564108 : j0lmn=jlmn*(jlmn-1)/2
514 1559602 : do ii=1,cplex
515 : gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
516 : & gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat) + &
517 : & enl_ptr2(j0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat)) * &
518 1516716 : & gx(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat)
519 : end do
520 : end do
521 : end do
522 : end do
523 : end do
524 :
525 : #ifdef HAVE_OPENMP_OFFLOAD
526 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
527 : !$OMP& MAP(to:enl_ptr2,atindx1,gx,gxfac_) &
528 : !$OMP& PRIVATE(idat,ispinor,ia,ispinor_index,index_enl) &
529 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
530 : #endif
531 47206 : do idat=1,ndat
532 79700 : do ispinor=1,nspinor
533 107874 : do ia=1,nincat
534 42886 : ispinor_index=ispinor+shift
535 42886 : index_enl=atindx1(iatm+ia)
536 : !$OMP PARALLEL DO PRIVATE(j0lmn,jlmn,ilmn,i0lmn,ii)
537 639488 : do jlmn=1,nlmn
538 564108 : j0lmn=jlmn*(jlmn-1)/2
539 4527666 : do ilmn=1,jlmn-1
540 11094594 : do ii=1,cplex
541 : gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat) + &
542 10530486 : & enl_ptr2(j0lmn+ilmn,index_enl,ispinor_index,min(ndat_enl,idat)) * gx(ii,ibeg+ilmn+(ia-1)*nlmn,ispinor,idat)
543 : end do
544 : end do
545 606994 : if(jlmn<nlmn) then
546 4484780 : do ilmn=jlmn+1,nlmn
547 3963558 : i0lmn=(ilmn*(ilmn-1)/2)
548 11051708 : do ii=1,cplex
549 : gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(ii,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat) &
550 : & + enl_ptr2(i0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat)) &
551 10530486 : & * gx(ii,ibeg+ilmn+(ia-1)*nlmn,ispinor,idat)
552 : end do
553 : end do
554 : end if
555 : end do
556 : end do
557 : end do
558 : end do
559 : endif
560 :
561 :
562 : ! 2-Enl is complex ===== D^ss'_ij=D^s's_ji^*
563 : else
564 0 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
565 :
566 0 : if (nspinortot==1) then ! -------------> NO SPINORS
567 0 : if(paw_opt==2) then
568 : #ifdef HAVE_OPENMP_OFFLOAD
569 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
570 : !$OMP& MAP(to:gxfac_,gx,gxi,atindx1,gxj,sij,enl_ptr2,lambda) &
571 : !$OMP& PRIVATE(idat,ia) &
572 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
573 : #endif
574 0 : do idat=1,ndat
575 0 : do ia=1,nincat
576 : !$OMP PARALLEL DO PRIVATE(index_enl,jlmn,j0lmn,enl_,gxj,ilmn,i0lmn,gxi)
577 0 : do jlmn=1,nlmn
578 0 : index_enl=atindx1(iatm+ia)
579 0 : j0lmn=jlmn*(jlmn-1)/2
580 0 : enl_(1)=enl_ptr2(2*j0lmn+jlmn-1,index_enl,1,min(ndat_enl,idat))-lambda(idat)*sij(j0lmn+jlmn)
581 0 : gxj(1 )=gx(1 ,jlmn+(ia-1)*nlmn+ibeg,1,idat)
582 0 : gxj(cplex)=gx(cplex,jlmn+(ia-1)*nlmn+ibeg,1,idat)
583 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
584 0 : & gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(1)
585 0 : if (cplex==2) then
586 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
587 0 : & gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(2)
588 : end if
589 0 : do ilmn=1,jlmn-1
590 0 : enl_(1)=enl_ptr2(2*j0lmn+ilmn-1,index_enl,1,min(ndat_enl,idat))
591 0 : enl_(2)=enl_ptr2(2*j0lmn+ilmn ,index_enl,1,min(ndat_enl,idat))
592 0 : enl_(1)=enl_(1)-lambda(idat)*sij(j0lmn+ilmn)
593 0 : gxi(1 )=gx(1 ,ilmn+(ia-1)*nlmn+ibeg,1,idat)
594 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,1,idat)
595 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
596 0 : & gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(1)
597 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
598 0 : & gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxi(1)
599 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat) = &
600 0 : & gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(1)
601 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat) = &
602 0 : & gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxj(1)
603 0 : if (cplex==2) then
604 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
605 0 : & gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxi(2)
606 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
607 0 : & gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(2)
608 : gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat) = &
609 0 : & gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxj(2)
610 : gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat) = &
611 0 : & gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(2)
612 : end if
613 : end do
614 0 : if(jlmn<nlmn) then
615 0 : do ilmn=jlmn+1,nlmn
616 0 : i0lmn=ilmn*(ilmn-1)/2
617 0 : enl_(1)=enl_ptr2(2*i0lmn+jlmn-1,index_enl,1,min(ndat_enl,idat))
618 0 : enl_(2)=enl_ptr2(2*i0lmn+jlmn ,index_enl,1,min(ndat_enl,idat))
619 0 : enl_(1)=enl_(1)-lambda(idat)*sij(i0lmn+jlmn)
620 0 : gxi(1 )=gx(1 ,ilmn+(ia-1)*nlmn+ibeg,1,idat)
621 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,1,idat)
622 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
623 0 : & gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(1)
624 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)= &
625 0 : & gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxi(1)
626 0 : if (cplex==2) then
627 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
628 0 : & gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxi(2)
629 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat) = &
630 0 : & gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(2)
631 : end if
632 : end do
633 : end if
634 : end do
635 : end do
636 : end do
637 : else
638 : #ifdef HAVE_OPENMP_OFFLOAD
639 : !$OMP TARGET TEAMS DISTRIBUTE &
640 : !$OMP& MAP(to:gxfac_,gx,gxi,atindx1,gxj,enl_ptr2) &
641 : !$OMP& PRIVATE(idat) &
642 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
643 : #endif
644 0 : do idat=1,ndat
645 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,index_enl,jlmn,j0lmn,enl_,gxj,ilmn,i0lmn,jjlmn,ijlmn,gxi)
646 0 : do ia=1,nincat
647 0 : do jlmn=1,nlmn
648 0 : index_enl=atindx1(iatm+ia)
649 0 : j0lmn=jlmn*(jlmn-1)/2
650 0 : jjlmn=j0lmn+jlmn
651 0 : enl_(1)=enl_ptr2(2*jjlmn-1,index_enl,1,min(ndat_enl,idat))
652 0 : gxj(1 )=gx(1 ,jlmn+(ia-1)*nlmn+ibeg,1,idat)
653 0 : gxj(cplex)=gx(cplex,jlmn+(ia-1)*nlmn+ibeg,1,idat)
654 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(1)
655 0 : if (cplex==2) gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(2)
656 0 : do ilmn=1,jlmn-1
657 0 : ijlmn=j0lmn+ilmn
658 0 : enl_(1)=enl_ptr2(2*ijlmn-1,index_enl,1,min(ndat_enl,idat))
659 0 : enl_(2)=enl_ptr2(2*ijlmn ,index_enl,1,min(ndat_enl,idat))
660 0 : gxi(1 )=gx(1 ,ilmn+(ia-1)*nlmn+ibeg,1,idat)
661 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,1,idat)
662 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(1)
663 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxi(1)
664 0 : if (cplex==2) then
665 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxi(2)
666 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(2)
667 : end if
668 : end do
669 0 : if(jlmn<nlmn) then
670 0 : do ilmn=jlmn+1,nlmn
671 0 : i0lmn=ilmn*(ilmn-1)/2
672 0 : ijlmn=i0lmn+jlmn
673 0 : enl_(1)=enl_ptr2(2*ijlmn-1,index_enl,1,min(ndat_enl,idat))
674 0 : enl_(2)=enl_ptr2(2*ijlmn ,index_enl,1,min(ndat_enl,idat))
675 0 : gxi(1 )=gx(1 ,ilmn+(ia-1)*nlmn+ibeg,1,idat)
676 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,1,idat)
677 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(1)
678 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxi(1)
679 0 : if (cplex==2) then
680 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxi(2)
681 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxi(2)
682 : end if
683 : end do
684 : end if
685 : end do
686 : end do
687 : end do
688 : end if
689 :
690 : else ! -------------> SPINORIAL CASE
691 :
692 : ! === Diagonal term(s) (up-up, down-down)
693 :
694 : #ifdef HAVE_OPENMP_OFFLOAD
695 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
696 : !$OMP& MAP(to:gxfac_,gx,enl_ptr,atindx1,sij) PRIVATE(idat,ispinor) &
697 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
698 : #endif
699 0 : do idat=1,ndat
700 0 : do ispinor=1,nspinor
701 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,ispinor_index,index_enl,j0lmn,jjlmn,i0lmn,ijlmn,gxi,gxj,enl_)
702 0 : do ia=1,nincat
703 0 : do jlmn=1,nlmn
704 0 : ispinor_index=ispinor+shift
705 0 : index_enl=atindx1(iatm+ia)
706 0 : j0lmn=jlmn*(jlmn-1)/2
707 0 : jjlmn=j0lmn+jlmn
708 0 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
709 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(jjlmn)
710 0 : gxj(1) =gx(1 ,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
711 0 : gxj(cplex)=gx(cplex,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
712 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxj(1)
713 0 : if (cplex==2) then
714 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxj(2)
715 : end if
716 0 : do ilmn=1,jlmn-1
717 0 : ijlmn=j0lmn+ilmn
718 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,ispinor_index)
719 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,ispinor_index)
720 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(ijlmn)
721 0 : gxi(1) =gx(1 ,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
722 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
723 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxi(1)
724 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)-enl_(2)*gxi(1)
725 0 : if (cplex==2) then
726 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(2)*gxi(2)
727 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxi(2)
728 : end if
729 : end do
730 : end do
731 : end do
732 : end do
733 : end do
734 : #ifdef HAVE_OPENMP_OFFLOAD
735 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
736 : !$OMP& MAP(to:gxfac_,gx,enl_ptr,atindx1) PRIVATE(idat,ispinor) &
737 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
738 : #endif
739 0 : do idat=1,ndat
740 0 : do ispinor=1,nspinor
741 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,ilmn,ispinor_index,index_enl,i0lmn,ijlmn,gxi,gxj,enl_)
742 0 : do ia=1,nincat
743 0 : do jlmn=1,nlmn-1
744 0 : do ilmn=jlmn+1,nlmn
745 0 : ispinor_index=ispinor+shift
746 0 : index_enl=atindx1(iatm+ia)
747 0 : i0lmn=ilmn*(ilmn-1)/2
748 0 : ijlmn=i0lmn+jlmn
749 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,ispinor_index)
750 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,ispinor_index)
751 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(ijlmn)
752 0 : gxi(1) =gx(1 ,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
753 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
754 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxi(1)
755 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(2)*gxi(1)
756 0 : if (cplex==2) then
757 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)-enl_(2)*gxi(2)
758 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxi(2)
759 : end if
760 : end do
761 : end do
762 : end do
763 : end do
764 : end do
765 : end if !nspinortot
766 : end if !complex_enl
767 :
768 : ! === Off-diagonal term(s) (up-down, down-up)
769 :
770 : ! --- No parallelization over spinors ---
771 15296 : if (nspinortot==2.and.nspinor==nspinortot) then
772 0 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
773 0 : ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex)!")
774 : #ifdef HAVE_OPENMP_OFFLOAD
775 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
776 : !$OMP& MAP(to:gxfac_,gx,enl_ptr,atindx1) PRIVATE(idat,ispinor) &
777 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
778 : #endif
779 0 : do idat=1,ndat
780 0 : do ispinor=1,nspinortot
781 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,jspinor,index_enl,j0lmn,jjlmn,i0lmn,ijlmn,gxi,enl_)
782 0 : do ia=1,nincat
783 0 : do jlmn=1,nlmn
784 0 : jspinor=3-ispinor
785 0 : index_enl=atindx1(iatm+ia)
786 0 : j0lmn=jlmn*(jlmn-1)/2
787 0 : jjlmn=j0lmn+jlmn
788 0 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,2+ispinor )
789 0 : enl_(2)=enl_ptr(2*jjlmn ,index_enl,2+ispinor )
790 0 : gxi(1) =gx(1 ,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
791 0 : gxi(cplex)=gx(cplex,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
792 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxi(1)
793 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)-enl_(2)*gxi(1)
794 0 : if (cplex==2) then
795 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(2)*gxi(2)
796 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxi(2)
797 : end if
798 0 : do ilmn=1,jlmn-1
799 0 : j0lmn=jlmn*(jlmn-1)/2
800 0 : ijlmn=j0lmn+ilmn
801 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,2+ispinor)
802 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,2+ispinor)
803 0 : gxi(1) =gx(1 ,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
804 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
805 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxi(1)
806 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)-enl_(2)*gxi(1)
807 0 : if (cplex==2) then
808 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(2)*gxi(2)
809 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxi(2)
810 : end if
811 : end do
812 : end do
813 : end do
814 : end do
815 : end do
816 : #ifdef HAVE_OPENMP_OFFLOAD
817 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
818 : !$OMP& MAP(to:gxfac_,gx,enl_ptr,atindx1) PRIVATE(idat,ispinor) &
819 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
820 : #endif
821 0 : do idat=1,ndat
822 0 : do ispinor=1,nspinortot
823 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,ilmn,jspinor,index_enl,i0lmn,ijlmn,gxi,enl_)
824 0 : do ia=1,nincat
825 0 : do jlmn=1,nlmn-1
826 0 : do ilmn=jlmn+1,nlmn
827 0 : jspinor=3-ispinor
828 0 : index_enl=atindx1(iatm+ia)
829 0 : i0lmn=ilmn*(ilmn-1)/2
830 0 : ijlmn=i0lmn+jlmn
831 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,2+ispinor)
832 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,2+ispinor)
833 0 : gxi(1) =gx(1 ,ilmn+(ia-1)*nlmn+ibeg,jspinor,idat)
834 0 : gxi(cplex)=gx(cplex,ilmn+(ia-1)*nlmn+ibeg,jspinor,idat)
835 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxi(1)
836 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(2)*gxi(1)
837 0 : if (cplex==2) then
838 0 : gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(1,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)-enl_(2)*gxi(2)
839 0 : gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=gxfac_(2,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxi(2)
840 : end if
841 : end do
842 : end do
843 : end do
844 : end do
845 : end do
846 :
847 : ! --- Parallelization over spinors ---
848 15296 : else if (nspinortot==2.and.nspinor/=nspinortot) then
849 0 : ABI_BUG("npspinor==2 not supported with OpenMP GPU")
850 : end if
851 :
852 : end if !paw_opt
853 :
854 :
855 : !Accumulate dgxdtfac related to nonlocal operator (Norm-conserving)
856 : !-------------------------------------------------------------------
857 18496 : if (optder>=1.and.paw_opt==0) then
858 : !Enl is E(Kleinman-Bylander)
859 0 : ABI_CHECK(cplex_enl==1,"BUG: invalid cplex_enl/=1!")
860 0 : ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!")
861 :
862 0 : if (lmax_so == 0) then
863 :
864 : ! NC+SR ---
865 : #ifdef HAVE_OPENMP_OFFLOAD
866 : !$OMP TARGET TEAMS DISTRIBUTE &
867 : !$OMP& MAP(to:dgxdtfac_,dgxdt,enl_ptr2,indlmn) &
868 : !$OMP& PRIVATE(idat,ispinor,ispinor_index) &
869 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
870 : #endif
871 0 : do idat=1,ndat
872 0 : do ispinor=1,nspinor
873 0 : ispinor_index = ispinor + shift
874 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(ia,ilmn,mu,ii)
875 0 : do ia=1,nincat
876 0 : do ilmn=1,nlmn
877 0 : do mu=1,ndgxdtfac
878 0 : do ii=1,cplex
879 0 : if (indlmn(6,ilmn)==2) then
880 0 : dgxdtfac_(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat) = zero
881 : else
882 : dgxdtfac_(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
883 : & enl_ptr2(indlmn(5,ilmn),itypat,ispinor_index,min(ndat_enl,idat)) &
884 0 : & *dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
885 : end if
886 : end do
887 : end do
888 : end do
889 : end do
890 : end do
891 : end do
892 :
893 : else
894 :
895 : ! NC+SR+SO: real-Ylm L.S coupling ---
896 : #ifdef HAVE_OPENMP_OFFLOAD
897 : !$OMP TARGET TEAMS DISTRIBUTE &
898 : !$OMP& MAP(to:dgxdtfac_,dgxdt,enl_ptr2,indlmn,ls_ylm_so_) &
899 : !$OMP& PRIVATE(idat) &
900 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
901 : #endif
902 0 : do idat=1,ndat
903 : !$OMP PARALLEL DO COLLAPSE(2) &
904 : !$OMP& PRIVATE(ia,ilmn,iln,ekb_so_1,ekb_so_2,ll_so,ilm,jlmn,jlm,klm_so,sign_so,ls_uu_im,ls_ud_re,ls_ud_im,mu)
905 0 : do ia=1,nincat
906 0 : do ilmn=1,nlmn
907 0 : iln = indlmn(5,ilmn)
908 :
909 0 : if (indlmn(6,ilmn)/=2) then
910 0 : ekb_so_1 = enl_ptr2(iln,itypat,1,min(ndat_enl,idat))
911 0 : ekb_so_2 = enl_ptr2(iln,itypat,2,min(ndat_enl,idat))
912 0 : do mu=1,ndgxdtfac
913 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)= &
914 0 : & ekb_so_1*dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
915 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)= &
916 0 : & ekb_so_1*dgxdt(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
917 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)= &
918 0 : & ekb_so_2*dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)
919 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)= &
920 0 : & ekb_so_2*dgxdt(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)
921 : end do
922 :
923 : else
924 :
925 0 : do mu=1,ndgxdtfac
926 0 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat) = zero
927 0 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat) = zero
928 0 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat) = zero
929 0 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat) = zero
930 : end do
931 :
932 0 : ekb_so_1 = enl_ptr2(iln,itypat,1,min(ndat_enl,idat))
933 0 : if (abs(ekb_so_1)<tol16) cycle
934 0 : ll_so = indlmn(1,ilmn)
935 0 : ilm = indlmn(4,ilmn)
936 0 : do jlmn=1,nlmn
937 0 : if (indlmn(6,jlmn)/=2) cycle
938 0 : if (indlmn(1,jlmn)/=ll_so) cycle
939 0 : if (indlmn(3,jlmn)/=indlmn(3,ilmn)) cycle
940 0 : jlm = indlmn(4,jlmn)
941 0 : if (ilm<=jlm) then
942 0 : klm_so = jlm*(jlm-1)/2 + ilm
943 0 : sign_so = 1
944 : else
945 0 : klm_so = ilm*(ilm-1)/2 + jlm
946 0 : sign_so = -1
947 : end if
948 0 : ls_uu_im = sign_so * ls_ylm_so_(2,klm_so,1)
949 0 : ls_ud_re = sign_so * ls_ylm_so_(1,klm_so,2)
950 0 : ls_ud_im = sign_so * ls_ylm_so_(2,klm_so,2)
951 :
952 0 : do mu=1,ndgxdtfac
953 : ! up-up
954 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
955 0 : & - ekb_so_1*ls_uu_im*dgxdt(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)
956 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
957 0 : & + ekb_so_1*ls_uu_im*dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)
958 : ! up-dn
959 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
960 0 : & + ekb_so_1*(ls_ud_re*dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,2,idat) - ls_ud_im*dgxdt(2,mu,jlmn+(ia-1)*nlmn+ibeg,2,idat))
961 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat) &
962 0 : & + ekb_so_1*(ls_ud_re*dgxdt(2,mu,jlmn+(ia-1)*nlmn+ibeg,2,idat) + ls_ud_im*dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,2,idat))
963 : ! dn-up
964 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)=dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
965 0 : & + ekb_so_1*(-ls_ud_re*dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat) - ls_ud_im*dgxdt(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat))
966 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)=dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
967 0 : & + ekb_so_1*(-ls_ud_re*dgxdt(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat) + ls_ud_im*dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat))
968 : ! dn-dn
969 : dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)=dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
970 0 : & + ekb_so_1*ls_uu_im*dgxdt(2,mu,jlmn+(ia-1)*nlmn+ibeg,2,idat)
971 : dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat)=dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,2,idat) &
972 0 : & - ekb_so_1*ls_uu_im*dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,2,idat)
973 : end do ! mu
974 : end do ! jlmn
975 : end if ! indlmn(:,6)==2
976 : end do ! ilmn
977 : end do ! ia
978 : !$OMP END PARALLEL DO
979 :
980 : end do ! idat
981 : end if ! NC+SO
982 : end if ! NC
983 :
984 18496 : if (lmax_so > 0) then
985 : #ifdef HAVE_OPENMP_OFFLOAD
986 : !$OMP TARGET EXIT DATA MAP(delete:ls_ylm_so_data) IF(gpu_option==ABI_GPU_OPENMP)
987 : #endif
988 0 : ABI_FREE(ls_ylm_so_data)
989 : end if
990 :
991 : !Accumulate dgxdtfac related to nonlocal operator (PAW)
992 : !-------------------------------------------------------------------
993 18496 : if (optder>=1.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then
994 : !Enl is psp strength Dij or (Dij-lambda.Sij)
995 :
996 : ! === Diagonal term(s) (up-up, down-down)
997 :
998 : ! 1-Enl is real
999 0 : if (cplex_enl==1) then
1000 0 : if (paw_opt/=2) then
1001 : #ifdef HAVE_OPENMP_OFFLOAD
1002 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1003 : !$OMP& MAP(to:dgxdtfac_,enl_ptr2,atindx1,dgxdt) &
1004 : !$OMP& PRIVATE(idat,ispinor,ia) &
1005 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1006 : #endif
1007 0 : do idat=1,ndat
1008 0 : do ispinor=1,nspinor
1009 0 : do ia=1,nincat
1010 : !$OMP PARALLEL DO &
1011 : !$OMP& PRIVATE(ispinor_index,index_enl,j0lmn,i0lmn,jlmn,ilmn,mu,ii)
1012 0 : do jlmn=1,nlmn
1013 0 : ispinor_index=ispinor+shift
1014 0 : index_enl=atindx1(iatm+ia)
1015 0 : j0lmn=jlmn*(jlmn-1)/2
1016 0 : do mu=1,ndgxdtfac
1017 0 : do ii=1,cplex
1018 : dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1019 : & dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1020 : & + enl_ptr2(j0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat)) &
1021 0 : & * dgxdt(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1022 : end do
1023 : end do
1024 0 : do ilmn=1,jlmn-1
1025 0 : do mu=1,ndgxdtfac
1026 0 : do ii=1,cplex
1027 : dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1028 : & dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1029 : & + enl_ptr2(j0lmn+ilmn,index_enl,ispinor_index,min(ndat_enl,idat)) &
1030 0 : & * dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1031 : end do
1032 : end do
1033 : end do
1034 0 : if(jlmn<nlmn) then
1035 0 : do ilmn=jlmn+1,nlmn
1036 0 : do mu=1,ndgxdtfac
1037 0 : do ii=1,cplex
1038 0 : i0lmn=ilmn*(ilmn-1)/2
1039 : dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1040 : & dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1041 : & + enl_ptr2(i0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat)) &
1042 0 : & * dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1043 : end do
1044 : end do
1045 : end do
1046 : end if
1047 : end do
1048 : end do
1049 : end do
1050 : end do
1051 : else
1052 : #ifdef HAVE_OPENMP_OFFLOAD
1053 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1054 : !$OMP& MAP(to:dgxdtfac_,enl_ptr2,atindx1,dgxdt,sij,lambda) &
1055 : !$OMP& PRIVATE(idat,ispinor) &
1056 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1057 : #endif
1058 0 : do idat=1,ndat
1059 0 : do ispinor=1,nspinor
1060 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ispinor_index,ia,index_enl,jlmn,j0lmn,ilmn,i0lmn,ii)
1061 0 : do ia=1,nincat
1062 0 : do jlmn=1,nlmn
1063 0 : ispinor_index=ispinor+shift
1064 0 : index_enl=atindx1(iatm+ia)
1065 0 : j0lmn=jlmn*(jlmn-1)/2
1066 0 : do mu=1,ndgxdtfac
1067 0 : do ii=1,cplex
1068 : dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1069 : & dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1070 : & + (enl_ptr2(j0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat))-lambda(idat)*sij(j0lmn+jlmn)) &
1071 0 : & * dgxdt(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1072 : end do
1073 : end do
1074 0 : do ilmn=1,jlmn-1
1075 0 : do mu=1,ndgxdtfac
1076 0 : do ii=1,cplex
1077 : dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1078 : & dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1079 : & + (enl_ptr2(j0lmn+ilmn,index_enl,ispinor_index,min(ndat_enl,idat))-lambda(idat)*sij(j0lmn+ilmn)) &
1080 0 : & * dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1081 : end do
1082 : end do
1083 : end do
1084 0 : if(jlmn<nlmn) then
1085 0 : do ilmn=jlmn+1,nlmn
1086 0 : i0lmn=ilmn*(ilmn-1)/2
1087 0 : do mu=1,ndgxdtfac
1088 0 : do ii=1,cplex
1089 : dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1090 : & dgxdtfac_(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1091 : & + (enl_ptr2(i0lmn+jlmn,index_enl,ispinor_index,min(ndat_enl,idat))-lambda(idat)*sij(i0lmn+jlmn)) &
1092 0 : & * dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1093 : end do
1094 : end do
1095 : end do
1096 : end if
1097 : end do
1098 : end do
1099 : end do
1100 : end do
1101 : end if
1102 :
1103 : ! 2-Enl is complex ===== D^ss'_ij=D^s's_ji^*
1104 : else
1105 0 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
1106 :
1107 0 : if (nspinortot==1) then ! -------------> NO SPINORS
1108 :
1109 : #ifdef HAVE_OPENMP_OFFLOAD
1110 : !$OMP TARGET TEAMS DISTRIBUTE &
1111 : !$OMP& MAP(to:dgxdtfac_,enl_,atindx1,dgxdt,sij,lambda,enl_ptr2,gxfi,gxj) &
1112 : !$OMP& PRIVATE(idat) &
1113 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1114 : #endif
1115 0 : do idat=1,ndat
1116 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,ilmn,i0lmn,ijlmn,gxfi,gxj,mu,cplex_,enl_)
1117 0 : do ia=1,nincat
1118 0 : do jlmn=1,nlmn
1119 0 : index_enl=atindx1(iatm+ia)
1120 0 : j0lmn=jlmn*(jlmn-1)/2
1121 0 : jjlmn=j0lmn+jlmn
1122 0 : enl_(1)=enl_ptr2(2*jjlmn-1,index_enl,1,min(ndat_enl,idat))
1123 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(jjlmn)
1124 0 : do mu=1,ndgxdtfac
1125 0 : if(cplex_dgxdt(mu)==2)then
1126 0 : cplex_ = 2 ; gxj(1) = zero ; gxj(2) = dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)
1127 : else
1128 0 : cplex_ = cplex ;
1129 0 : gxj(1 )=dgxdt(1 ,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)
1130 0 : gxj(cplex)=dgxdt(cplex,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)
1131 : end if
1132 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(1)
1133 0 : if (cplex_==2) then
1134 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxj(2)
1135 : end if
1136 : end do
1137 0 : do ilmn=1,jlmn-1
1138 0 : ijlmn=j0lmn+ilmn
1139 0 : enl_(1)=enl_ptr2(2*ijlmn-1,index_enl,1,min(ndat_enl,idat))
1140 0 : enl_(2)=enl_ptr2(2*ijlmn ,index_enl,1,min(ndat_enl,idat))
1141 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(ijlmn)
1142 0 : do mu=1,ndgxdtfac
1143 0 : if(cplex_dgxdt(mu)==2)then
1144 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
1145 : else
1146 0 : cplex_ = cplex ;
1147 0 : gxfi(1 )=dgxdt(1 ,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
1148 0 : gxfi(cplex)=dgxdt(cplex,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
1149 : end if
1150 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxfi(1)
1151 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxfi(1)
1152 0 : if (cplex_==2) then
1153 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxfi(2)
1154 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxfi(2)
1155 : end if
1156 : end do
1157 : end do
1158 0 : if(jlmn<nlmn) then
1159 0 : do ilmn=jlmn+1,nlmn
1160 0 : i0lmn=ilmn*(ilmn-1)/2
1161 0 : ijlmn=i0lmn+jlmn
1162 0 : enl_(1)=enl_ptr2(2*ijlmn-1,index_enl,1,min(ndat_enl,idat))
1163 0 : enl_(2)=enl_ptr2(2*ijlmn ,index_enl,1,min(ndat_enl,idat))
1164 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(ijlmn)
1165 0 : do mu=1,ndgxdtfac
1166 0 : if(cplex_dgxdt(mu)==2)then
1167 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
1168 : else
1169 0 : cplex_ = cplex ;
1170 0 : gxfi(1 )=dgxdt(1 ,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
1171 0 : gxfi(cplex)=dgxdt(cplex,mu,ilmn+(ia-1)*nlmn+ibeg,1,idat)
1172 : end if
1173 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxfi(1)
1174 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(2)*gxfi(1)
1175 0 : if (cplex_==2) then
1176 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)-enl_(2)*gxfi(2)
1177 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,1,idat)+enl_(1)*gxfi(2)
1178 : end if
1179 : end do
1180 : end do
1181 : end if
1182 : end do
1183 : end do
1184 : end do
1185 : else ! -------------> SPINORIAL CASE
1186 :
1187 : ! === Diagonal term(s) (up-up, down-down)
1188 :
1189 : #ifdef HAVE_OPENMP_OFFLOAD
1190 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1191 : !$OMP& MAP(to:atindx1,dgxdtfac_,enl_ptr,dgxdt,sij,lambda) PRIVATE(idat,ispinor) &
1192 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1193 : #endif
1194 0 : do idat=1,ndat
1195 0 : do ispinor=1,nspinor
1196 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,ispinor_index,index_enl,j0lmn,jjlmn,ilmn,ijlmn,cplex_,enl_,gxj,gxfi)
1197 0 : do ia=1,nincat
1198 0 : do jlmn=1,nlmn
1199 0 : ispinor_index = ispinor + shift
1200 0 : index_enl=atindx1(iatm+ia)
1201 0 : j0lmn=jlmn*(jlmn-1)/2
1202 0 : jjlmn=j0lmn+jlmn
1203 0 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
1204 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(jjlmn)
1205 0 : do mu=1,ndgxdtfac
1206 0 : if(cplex_dgxdt(mu)==2)then
1207 0 : cplex_ = 2
1208 0 : gxj(1) = zero ; gxj(2) = dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1209 : else
1210 0 : cplex_ = cplex
1211 0 : gxj(1 )=dgxdt(1 ,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1212 0 : gxj(cplex)=dgxdt(cplex,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1213 : end if
1214 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxj(1)
1215 0 : if (cplex_==2) then
1216 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxj(2)
1217 : end if
1218 : end do
1219 0 : do ilmn=1,jlmn-1
1220 0 : ijlmn=j0lmn+ilmn
1221 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,ispinor_index)
1222 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,ispinor_index)
1223 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(ijlmn)
1224 0 : do mu=1,ndgxdtfac
1225 0 : if(cplex_dgxdt(mu)==2)then
1226 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1227 : else
1228 0 : cplex_ = cplex
1229 0 : gxfi(1) =dgxdt(1 ,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1230 0 : gxfi(cplex)=dgxdt(cplex,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1231 : end if
1232 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxfi(1)
1233 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)-enl_(2)*gxfi(1)
1234 0 : if (cplex_==2) then
1235 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(2)*gxfi(2)
1236 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxfi(2)
1237 : end if
1238 : end do
1239 : end do
1240 : end do
1241 : end do
1242 : end do
1243 : end do
1244 : #ifdef HAVE_OPENMP_OFFLOAD
1245 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1246 : !$OMP& MAP(to:atindx1,dgxdtfac_,enl_ptr,dgxdt,sij,lambda) PRIVATE(idat,ispinor) &
1247 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1248 : #endif
1249 0 : do idat=1,ndat
1250 0 : do ispinor=1,nspinor
1251 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,ispinor_index,index_enl,i0lmn,ilmn,ijlmn,cplex_,enl_,gxfi)
1252 0 : do ia=1,nincat
1253 0 : do jlmn=1,nlmn-1
1254 0 : do ilmn=jlmn+1,nlmn
1255 0 : ispinor_index = ispinor + shift
1256 0 : index_enl=atindx1(iatm+ia)
1257 0 : i0lmn=ilmn*(ilmn-1)/2
1258 0 : ijlmn=i0lmn+jlmn
1259 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,ispinor_index)
1260 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,ispinor_index)
1261 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda(idat)*sij(ijlmn)
1262 0 : do mu=1,ndgxdtfac
1263 0 : if(cplex_dgxdt(mu)==2)then
1264 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1265 : else
1266 0 : cplex_ = cplex ;
1267 0 : gxfi(1 )=dgxdt(1 ,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1268 0 : gxfi(cplex)=dgxdt(cplex,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1269 : end if
1270 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxfi(1)
1271 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(2)*gxfi(1)
1272 0 : if (cplex_==2) then
1273 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)-enl_(2)*gxfi(2)
1274 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxfi(2)
1275 : end if
1276 : end do
1277 : end do
1278 : end do
1279 : end do
1280 : end do
1281 : end do
1282 : end if !nspinortot
1283 : end if !complex
1284 :
1285 : ! === Off-diagonal term(s) (up-down, down-up)
1286 :
1287 : ! --- No parallelization over spinors ---
1288 0 : if (nspinortot==2.and.nspinor==nspinortot) then
1289 0 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
1290 0 : ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
1291 : #ifdef HAVE_OPENMP_OFFLOAD
1292 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1293 : !$OMP& MAP(to:atindx1,dgxdtfac_,enl_ptr,dgxdt) PRIVATE(idat,ispinor) &
1294 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1295 : #endif
1296 0 : do idat=1,ndat
1297 0 : do ispinor=1,nspinor
1298 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,jspinor,index_enl,j0lmn,jjlmn,ilmn,ijlmn,cplex_,enl_,gxj,gxfi)
1299 0 : do ia=1,nincat
1300 0 : do jlmn=1,nlmn
1301 0 : jspinor=3-ispinor
1302 0 : index_enl=atindx1(iatm+ia)
1303 0 : j0lmn=jlmn*(jlmn-1)/2
1304 0 : jjlmn=j0lmn+jlmn
1305 0 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,2+ispinor)
1306 0 : enl_(2)=enl_ptr(2*jjlmn ,index_enl,2+ispinor)
1307 0 : do mu=1,ndgxdtfac
1308 0 : if(cplex_dgxdt(mu)==2)then
1309 0 : cplex_ = 2 ;
1310 0 : gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1311 : else
1312 0 : cplex_ = cplex ;
1313 0 : gxfi(1) =dgxdt(1 ,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1314 0 : gxfi(cplex)=dgxdt(cplex,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1315 : end if
1316 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxfi(1)
1317 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)-enl_(2)*gxfi(1)
1318 0 : if (cplex_==2) then
1319 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(2)*gxfi(2)
1320 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxfi(2)
1321 : end if
1322 : end do
1323 0 : do ilmn=1,jlmn-1
1324 0 : ijlmn=j0lmn+ilmn
1325 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,2+ispinor)
1326 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,2+ispinor)
1327 0 : do mu=1,ndgxdtfac
1328 0 : if(cplex_dgxdt(mu)==2)then
1329 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1330 : else
1331 0 : cplex_ = cplex
1332 0 : gxfi(1) =dgxdt(1 ,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1333 0 : gxfi(cplex)=dgxdt(cplex,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1334 : end if
1335 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxfi(1)
1336 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)-enl_(2)*gxfi(1)
1337 0 : if (cplex_==2) then
1338 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(2)*gxfi(2)
1339 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,jspinor,idat)+enl_(1)*gxfi(2)
1340 : end if
1341 : end do !mu
1342 : end do !ilmn
1343 : end do !jmln
1344 : end do !ia
1345 : end do !ispinor
1346 : end do !idat
1347 : #ifdef HAVE_OPENMP_OFFLOAD
1348 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1349 : !$OMP& MAP(to:atindx1,dgxdtfac_,enl_ptr,dgxdt) PRIVATE(idat,ispinor) &
1350 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1351 : #endif
1352 0 : do idat=1,ndat
1353 0 : do ispinor=1,nspinor
1354 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,jspinor,index_enl,i0lmn,ilmn,ijlmn,cplex_,enl_,gxfi)
1355 0 : do ia=1,nincat
1356 0 : do jlmn=1,nlmn
1357 0 : do ilmn=jlmn+1,nlmn
1358 0 : jspinor=3-ispinor
1359 0 : index_enl=atindx1(iatm+ia)
1360 0 : i0lmn=ilmn*(ilmn-1)/2
1361 0 : ijlmn=i0lmn+jlmn
1362 0 : enl_(1)=enl_ptr(2*ijlmn-1,index_enl,2+ispinor)
1363 0 : enl_(2)=enl_ptr(2*ijlmn ,index_enl,2+ispinor)
1364 0 : do mu=1,ndgxdtfac
1365 0 : if(cplex_dgxdt(mu)==2)then
1366 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn+(ia-1)*nlmn+ibeg,jspinor,idat)
1367 : else
1368 0 : cplex_ = cplex
1369 0 : gxfi(1) =dgxdt(1 ,mu,ilmn+(ia-1)*nlmn+ibeg,jspinor,idat)
1370 0 : gxfi(cplex)=dgxdt(cplex,mu,ilmn+(ia-1)*nlmn+ibeg,jspinor,idat)
1371 : end if
1372 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxfi(1)
1373 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(2)*gxfi(1)
1374 0 : if (cplex_==2) then
1375 0 : dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(1,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)-enl_(2)*gxfi(2)
1376 0 : dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=dgxdtfac_(2,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)+enl_(1)*gxfi(2)
1377 : end if
1378 : end do !mu
1379 : end do !ilmn
1380 : end do !jmln
1381 : end do !ia
1382 : end do !ispinor
1383 : end do !idat
1384 :
1385 : ! --- Parallelization over spinors ---
1386 0 : else if (nspinortot==2.and.nspinor/=nspinortot) then
1387 0 : ABI_BUG("nspinor==2 not supported with OpenMP GPU")
1388 : end if !nspinortot
1389 :
1390 : end if ! pawopt & optder
1391 :
1392 : !End of loop when a exp(-iqR) phase is present
1393 : !------------------------------------------- ------------------------
1394 :
1395 : !When iphase=1, gxfac and gxfac_ point to the same memory space
1396 : !When iphase=2, we add i.gxfac_ to gxfac
1397 36992 : if (iphase==2) then
1398 : #ifdef HAVE_OPENMP_OFFLOAD
1399 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1400 : !$OMP& PRIVATE(idat,ia) MAP(to:gxfac,gxfac_) &
1401 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1402 : #endif
1403 0 : do idat=1,ndat
1404 0 : do ispinor=1,nspinor
1405 0 : do ia=1,nincat
1406 : !$OMP PARALLEL DO PRIVATE(ilmn)
1407 0 : do ilmn=1,nlmn
1408 : gxfac(1,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1409 0 : & gxfac(1,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)-gxfac_(2,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1410 : gxfac(2,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1411 0 : & gxfac(2,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)+gxfac_(1,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1412 : end do
1413 : end do
1414 : end do
1415 : end do
1416 0 : if (optder>=1) then
1417 : #ifdef HAVE_OPENMP_OFFLOAD
1418 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1419 : !$OMP& PRIVATE(idat,ia,ilmn,mu) MAP(to:dgxdtfac,dgxdtfac_) &
1420 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1421 : #endif
1422 0 : do idat=1,ndat
1423 0 : do ispinor=1,nspinor
1424 0 : do ia=1,nincat
1425 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ilmn,mu)
1426 0 : do ilmn=1,nlmn
1427 0 : do mu=1,ndgxdtfac
1428 : dgxdtfac(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1429 0 : & dgxdtfac(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)-dgxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1430 : dgxdtfac(2,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1431 0 : & dgxdtfac(2,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)+dgxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1432 : end do
1433 : end do
1434 : end do
1435 : end do
1436 : end do
1437 : end if
1438 0 : if (optder>=2) then
1439 : #ifdef HAVE_OPENMP_OFFLOAD
1440 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1441 : !$OMP& PRIVATE(idat,ia) MAP(to:d2gxdtfac,d2gxdtfac_) &
1442 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1443 : #endif
1444 0 : do idat=1,ndat
1445 0 : do ispinor=1,nspinor
1446 0 : do ia=1,nincat
1447 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ilmn,mu)
1448 0 : do ilmn=1,nlmn
1449 0 : do mu=1,nd2gxdtfac
1450 : d2gxdtfac(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1451 0 : & d2gxdtfac(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)-d2gxdtfac_(2,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1452 : d2gxdtfac(2,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1453 0 : & d2gxdtfac(2,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)+d2gxdtfac_(1,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1454 : end do
1455 : end do
1456 : end do
1457 : end do
1458 : end do
1459 : end if
1460 0 : call destroy_work_arrays(gpu_option)
1461 : end if
1462 :
1463 : !End loop over real/imaginary part of the exp(-iqR) phase
1464 : end do
1465 :
1466 :
1467 : !Accumulate gxfac related to overlap (Sij) (PAW)
1468 : !------------------------------------------- ------------------------
1469 18496 : if (paw_opt==3.or.paw_opt==4) then ! Use Sij, overlap contribution
1470 : #ifdef HAVE_OPENMP_OFFLOAD
1471 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1472 : !$OMP& MAP(to:sij,gx,gxfac_sij) &
1473 : !$OMP& PRIVATE(idat,ispinor) &
1474 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1475 : #endif
1476 43782 : do idat=1,ndat
1477 73844 : do ispinor=1,nspinor
1478 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,j0lmn,jjlmn,ilmn,i0lmn,ijlmn,ii)
1479 98850 : do ia=1,nincat
1480 592576 : do jlmn=1,nlmn
1481 523788 : j0lmn=jlmn*(jlmn-1)/2
1482 523788 : jjlmn=j0lmn+jlmn
1483 1400508 : do ii=1,cplex
1484 : gxfac_sij(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat)= &
1485 : gxfac_sij(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat) &
1486 1400508 : + sij(jjlmn) * gx(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat)
1487 : end do
1488 4282866 : do ilmn=1,jlmn-1
1489 3759078 : ijlmn=j0lmn+ilmn
1490 10481226 : do ii=1,cplex
1491 : gxfac_sij(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat)= &
1492 : gxfac_sij(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat) &
1493 9957438 : + sij(ijlmn) * gx(ii,ibeg+ilmn+(ia-1)*nlmn,ispinor,idat)
1494 : end do
1495 : end do
1496 562514 : if(jlmn<nlmn) then
1497 4244140 : do ilmn=jlmn+1,nlmn
1498 3759078 : i0lmn=ilmn*(ilmn-1)/2
1499 3759078 : ijlmn=i0lmn+jlmn
1500 10442500 : do ii=1,cplex
1501 : gxfac_sij(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat)=&
1502 : gxfac_sij(ii,ibeg+jlmn+(ia-1)*nlmn,ispinor,idat) &
1503 9957438 : + sij(ijlmn) * gx(ii,ibeg+ilmn+(ia-1)*nlmn,ispinor,idat)
1504 : end do
1505 : end do
1506 : end if
1507 : end do
1508 : end do
1509 : end do
1510 : end do
1511 : end if
1512 :
1513 : !Accumulate dgxdtfac related to overlap (Sij) (PAW)
1514 : !-------------------------------------------------------------------
1515 18496 : if (optder>=1.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution
1516 : #ifdef HAVE_OPENMP_OFFLOAD
1517 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1518 : !$OMP& MAP(to:sij,dgxdt,dgxdtfac_sij) &
1519 : !$OMP& PRIVATE(idat,ispinor) &
1520 : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1521 : #endif
1522 0 : do idat=1,ndat
1523 0 : do ispinor=1,nspinor
1524 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,jlmn,j0lmn,jjlmn,ilmn,i0lmn,ijlmn,ii)
1525 0 : do ia=1,nincat
1526 0 : do jlmn=1,nlmn
1527 0 : j0lmn=jlmn*(jlmn-1)/2
1528 0 : jjlmn=j0lmn+jlmn
1529 0 : do mu=1,ndgxdtfac
1530 0 : do ii=1,cplex
1531 : dgxdtfac_sij(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1532 : & dgxdtfac_sij(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1533 0 : & + sij(jjlmn) * dgxdt(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1534 : end do
1535 : end do
1536 0 : do ilmn=1,jlmn-1
1537 0 : ijlmn=j0lmn+ilmn
1538 0 : do mu=1,ndgxdtfac
1539 0 : do ii=1,cplex
1540 : dgxdtfac_sij(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1541 : & dgxdtfac_sij(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1542 0 : & + sij(ijlmn) * dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1543 : end do
1544 : end do
1545 : end do
1546 0 : if(jlmn<nlmn) then
1547 0 : do ilmn=jlmn+1,nlmn
1548 0 : i0lmn=ilmn*(ilmn-1)/2
1549 0 : ijlmn=i0lmn+jlmn
1550 0 : do mu=1,ndgxdtfac
1551 0 : do ii=1,cplex
1552 : dgxdtfac_sij(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)=&
1553 : & dgxdtfac_sij(ii,mu,jlmn+(ia-1)*nlmn+ibeg,ispinor,idat)&
1554 0 : & + sij(ijlmn) * dgxdt(ii,mu,ilmn+(ia-1)*nlmn+ibeg,ispinor,idat)
1555 : end do
1556 : end do
1557 : end do
1558 : end if
1559 : end do
1560 : end do
1561 : end do
1562 : end do
1563 : end if
1564 :
1565 18496 : end subroutine opernlc_ylm_allwf
1566 : !!***
1567 :
1568 : end module m_opernlc_ylm_allwf
1569 : !!***
|
