Branch data Line data Source code
1 : : !!****m* ABINIT/m_opernlc_ylm
2 : : !! NAME
3 : : !! m_opernlc_ylm
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
22 : :
23 : : use defs_basis
24 : : use m_errors
25 : : use m_abicore
26 : : use m_xmpi
27 : : use m_xomp
28 : :
29 : : use defs_abitypes, only : MPI_type
30 : : implicit none
31 : :
32 : : private
33 : : !!***
34 : :
35 : : public :: opernlc_ylm
36 : : public :: ls_ylm
37 : : !!***
38 : :
39 : : contains
40 : : !!***
41 : :
42 : : !!****f* ABINIT/opernlc_ylm
43 : : !! NAME
44 : : !! opernlc_ylm
45 : : !!
46 : : !! FUNCTION
47 : : !! * Operate with the non-local part of the hamiltonian,
48 : : !! in order to reduce projected scalars
49 : : !! * Operate with the non-local projectors and the overlap matrix,
50 : : !! in order to reduce projected scalars
51 : : !!
52 : : !! INPUTS
53 : : !! atindx1(natom)=index table for atoms (gives the absolute index of
54 : : !! an atom from its rank in a block of atoms)
55 : : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
56 : : !! 2 if <p_lmn|c> scalars are complex
57 : : !! cplex_dgxdt(ndgxdt) = used only when cplex = 1
58 : : !! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
59 : : !! cplex_enl=1 if enl factors are real, 2 if they are complex
60 : : !! cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
61 : : !! dgxdt(cplex,ndgxdt,nlmn,nincat)=grads of projected scalars (only if optder>0)
62 : : !! dimenl1,dimenl2=dimensions of enl (see enl)
63 : : !! dimekbq=1 if enl factors do not contain a exp(-iqR) phase, 2 is they do
64 : : !! enl(cplex_enl*dimenl1,dimenl2,nspinortot**2,dimekbq)=
65 : : !! ->Norm conserving : ==== when paw_opt=0 ====
66 : : !! (Real) Kleinman-Bylander energies (hartree)
67 : : !! dimenl1=lmnmax - dimenl2=ntypat
68 : : !! dimekbq is 2 if Enl contains a exp(-iqR) phase, 1 otherwise
69 : : !! ->PAW : ==== when paw_opt=1, 2 or 4 ====
70 : : !! (Real or complex, hermitian) Dij coefs to connect projectors
71 : : !! dimenl1=cplex_enl*lmnmax*(lmnmax+1)/2 - dimenl2=natom
72 : : !! These are complex numbers if cplex_enl=2
73 : : !! enl(:,:,1) contains Dij^up-up
74 : : !! enl(:,:,2) contains Dij^dn-dn
75 : : !! enl(:,:,3) contains Dij^up-dn (only if nspinor=2)
76 : : !! enl(:,:,4) contains Dij^dn-up (only if nspinor=2)
77 : : !! dimekbq is 2 if Dij contains a exp(-iqR) phase, 1 otherwise
78 : : !! gx(cplex,nlmn,nincat*abs(enl_opt))= projected scalars
79 : : !! iatm=absolute rank of first atom of the current block of atoms
80 : : !! indlmn(6,nlmn)= array giving l,m,n,lm,ln,s for i=lmn
81 : : !! itypat=type of atoms
82 : : !! lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
83 : : !! mpi_enreg=information about MPI parallelization
84 : : !! natom=number of atoms in cell
85 : : !! ndgxdt=second dimension of dgxdt
86 : : !! ndgxdtfac=second dimension of dgxdtfac
87 : : !! nincat=number of atoms in the subset here treated
88 : : !! nlmn=number of (l,m,n) numbers for current type of atom
89 : : !! nspinor= number of spinorial components of the wavefunctions (on current proc)
90 : : !! nspinortot=total number of spinorial components of the wavefunctions
91 : : !! optder=0=only gxfac is computed, 1=both gxfac and dgxdtfac are computed
92 : : !! 2=gxfac, dgxdtfac and d2gxdtfac are computed
93 : : !! paw_opt= define the nonlocal operator concerned with:
94 : : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
95 : : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
96 : : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
97 : : !! paw_opt=3 : PAW overlap matrix (Sij)
98 : : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
99 : : !! sij(nlm*(nlmn+1)/2)=overlap matrix components (only if paw_opt=2, 3 or 4)
100 : : !!
101 : : !! OUTPUT
102 : : !! if (paw_opt=0, 1, 2 or 4)
103 : : !! gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
104 : : !! if (paw_opt=3 or 4)
105 : : !! gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
106 : : !! if (optder==1.and.paw_opt=0, 1, 2 or 4)
107 : : !! dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Vnl (NL operator)
108 : : !! if (optder==1.and.paw_opt=3 or 4)
109 : : !! dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor)= gradients of gxfac related to Sij (overlap)
110 : : !!
111 : : !! NOTES
112 : : !! This routine operates for one type of atom, and within this given type of atom,
113 : : !! for a subset of at most nincat atoms.
114 : : !!
115 : : !! About the non-local factors symmetry:
116 : : !! - The lower triangular part of the Dij matrix can be deduced from the upper one
117 : : !! with the following relation: D^s2s1_ji = (D^s1s2_ij)^*
118 : : !! where s1,s2 are spinor components
119 : : !! - The Dij factors can contain a exp(-iqR) phase
120 : : !! This phase does not have to be included in the symmetry rule
121 : : !! For that reason, we first apply the real part (cos(qR).D^s1s2_ij)
122 : : !! then, we apply the imaginary part (-sin(qR).D^s1s2_ij)
123 : : !!
124 : : !! SOURCE
125 : :
126 : 35485077 : subroutine opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,&
127 : 35485077 : & dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,dimekbq,enl,&
128 : 35485077 : & gx,gxfac,gxfac_sij,iatm,indlmn,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,&
129 : 35485077 : & nd2gxdt,nd2gxdtfac,nincat,nlmn,nspinor,nspinortot,optder,paw_opt,sij)
130 : :
131 : : !Arguments ------------------------------------
132 : : !scalars
133 : : integer,intent(in) :: cplex,cplex_enl,cplex_fac,dimenl1,dimenl2,dimekbq,iatm,itypat
134 : : integer,intent(in) :: natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,nincat,nspinor,nspinortot,optder,paw_opt
135 : : integer,intent(inout) :: nlmn
136 : : real(dp) :: lambda
137 : : type(MPI_type) , intent(in) :: mpi_enreg
138 : : !arrays
139 : : integer,intent(in) :: atindx1(natom),indlmn(6,nlmn),cplex_dgxdt(ndgxdt),cplex_d2gxdt(nd2gxdt)
140 : : real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)
141 : : real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)
142 : : real(dp),intent(in),target :: enl(dimenl1,dimenl2,nspinortot**2,dimekbq)
143 : : real(dp),intent(inout) :: gx(cplex,nlmn,nincat,nspinor)
144 : : real(dp),intent(in) :: sij(((paw_opt+1)/3)*nlmn*(nlmn+1)/2)
145 : : real(dp),intent(out),target :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)
146 : : real(dp),intent(out) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3))
147 : : real(dp),intent(out),target :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor)
148 : : real(dp),intent(out) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat,nspinor*(paw_opt/3))
149 : : real(dp),intent(out),target :: gxfac(cplex_fac,nlmn,nincat,nspinor)
150 : : real(dp),intent(out) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
151 : :
152 : : !Local variables-------------------------------
153 : : !Arrays
154 : : !scalars
155 : : integer :: cplex_,ia,ierr,ijlmn,ijspin,ilm,ilmn,i0lmn,iln,index_enl,iphase,ispinor,ispinor_index
156 : : integer :: j0lmn,jilmn,jispin,jjlmn,jlm,jlmn,jspinor,jspinor_index,mu,shift
157 : : integer :: ll_so, klm_so, lmax_so, nlmso, sign_so
158 : : real(dp) :: sijr
159 : : real(dp) :: ekb_so, ls_uu_im, ls_ud_re, ls_ud_im
160 : 35485077 : real(dp), allocatable :: ls_ylm_so(:,:,:)
161 : : !arrays
162 : 70970154 : real(dp) :: enl_(2),gxfi(2),gxi(cplex),gxj(cplex)
163 : 35485077 : real(dp),allocatable :: d2gxdtfac_offdiag(:,:,:,:,:),dgxdtfac_offdiag(:,:,:,:,:)
164 : 35485077 : real(dp),allocatable :: gxfac_offdiag(:,:,:,:),gxfj(:,:)
165 : 35485077 : real(dp),pointer :: d2gxdtfac_(:,:,:,:,:),dgxdtfac_(:,:,:,:,:),gxfac_(:,:,:,:)
166 : 35485077 : real(dp),pointer :: enl_ptr(:,:,:)
167 : :
168 : : ! *************************************************************************
169 : :
170 : : DBG_ENTER("COLL")
171 : :
172 : : !Parallelization over spinors treatment
173 : 60248 : shift=0;if (mpi_enreg%paral_spinor==1) shift=mpi_enreg%me_spinor
174 : :
175 : : !When Enl factors contain a exp(-iqR) phase:
176 : : ! - We loop over the real and imaginary parts
177 : : ! - We need an additional memory space
178 [ + + ]: 76298294 : do iphase=1,dimekbq
179 [ + + ]: 40813217 : if (paw_opt==3) cycle
180 [ + + ]: 38506025 : if (iphase==1) then
181 : 33177885 : gxfac_ => gxfac ; dgxdtfac_ => dgxdtfac ; d2gxdtfac_ => d2gxdtfac
182 : : else
183 [ - + ]: 5328140 : ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac==1 when dimekbq=2!")
184 [ - + + - : 31968840 : ABI_MALLOC(gxfac_,(cplex_fac,nlmn,nincat,nspinor))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
185 [ + + ]: 5328140 : if (optder>=1) then
186 [ - + + - : 775908 : ABI_MALLOC(dgxdtfac_,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + ]
187 : : end if
188 [ - + ]: 5328140 : if (optder>=2) then
189 [ # # # # : 0 : ABI_MALLOC(d2gxdtfac_,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
190 : : end if
191 : : end if
192 [ + + + + : 1718605737 : gxfac_=zero
+ + + + ]
193 [ + + + + : 272789529 : if (optder>=1) dgxdtfac_=zero
+ + + + +
+ + + ]
194 [ + + + + : 108823286 : if (optder>=2) d2gxdtfac_=zero
+ + + + +
+ + + ]
195 : 38506025 : enl_ptr => enl(:,:,:,iphase)
196 : :
197 : : !NC+SO: precompute L.S matrix once (reused by gxfac, dgxdtfac, d2gxdtfac blocks below)
198 : 38506025 : lmax_so = 0
199 [ + + + + : 38506025 : if (paw_opt==0.and.nspinortot==2.and.nspinor==nspinortot) then
+ - ]
200 [ + + + + : 137970 : if (any(indlmn(6,1:nlmn)==2)) then
+ + ]
201 [ + + ]: 170820 : do ilmn=1,nlmn
202 [ + + ]: 170820 : if (indlmn(6,ilmn)==2) lmax_so = max(lmax_so, indlmn(1,ilmn))
203 : : end do
204 [ + - ]: 6570 : if (lmax_so > 0) then
205 : 6570 : nlmso = (lmax_so+1)**2*((lmax_so+1)**2+1)/2
206 [ - + + - : 26280 : ABI_MALLOC(ls_ylm_so,(2,nlmso,2))
+ - + - -
+ - + -
+ ]
207 : 6570 : call ls_ylm(ls_ylm_so, lmax_so)
208 : : end if
209 : : end if
210 : : end if
211 : :
212 : : !Accumulate gxfac related to non-local operator (Norm-conserving)
213 : : !-------------------------------------------------------------------
214 [ + + ]: 38506025 : if (paw_opt==0) then
215 : : !Enl is E(Kleinman-Bylander)
216 [ - + ]: 9449245 : ABI_CHECK(cplex_enl/=2,"BUG: invalid cplex_enl=2!")
217 [ - + ]: 9449245 : ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!")
218 : : !$OMP PARALLEL &
219 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_)
220 : : #ifdef FC_NVHPC
221 : : !FIXME Compiler bug since v24.1~24.9
222 : : !$OMP DO COLLAPSE(2)
223 : : #else
224 : : !$OMP DO COLLAPSE(3)
225 : : #endif
226 [ + + ]: 18911630 : do ispinor=1,nspinor
227 [ + + ]: 30758928 : do ia=1,nincat
228 [ + + ]: 144341753 : do ilmn=1,nlmn
229 [ + + ]: 123032070 : if (indlmn(6,ilmn)==2) cycle ! NC+SO: SO projectors handled separately below
230 : 122887530 : ispinor_index=ispinor+shift
231 : 122887530 : iln=indlmn(5,ilmn)
232 : 122887530 : enl_(1)=enl_ptr(iln,itypat,ispinor_index)
233 [ + + ]: 380132032 : gxfac_(1:cplex,ilmn,ia,ispinor)=enl_(1)*gx(1:cplex,ilmn,ia,ispinor)
234 : : end do
235 : : end do
236 : : end do
237 : : !$OMP END DO
238 : : !$OMP END PARALLEL
239 : : end if
240 : :
241 : : ! NC+SO: real-Ylm L.S coupling ---
242 [ + + ]: 38506025 : if (lmax_so > 0) then
243 [ + + ]: 13140 : do ia=1,nincat
244 [ + + ]: 177390 : do ilmn=1,nlmn
245 [ + + ]: 164250 : if (indlmn(6,ilmn)/=2) cycle
246 : 72270 : iln = indlmn(5,ilmn)
247 : 72270 : ekb_so = enl_ptr(iln,itypat,1)
248 [ - + ]: 72270 : if (abs(ekb_so)<tol16) cycle
249 : 72270 : ll_so = indlmn(1,ilmn)
250 : 72270 : ilm = indlmn(4,ilmn)
251 [ + + ]: 1885590 : do jlmn=1,nlmn
252 [ + + ]: 1806750 : if (indlmn(6,jlmn)/=2) cycle
253 [ + + ]: 794970 : if (indlmn(1,jlmn)/=ll_so) cycle
254 [ + + ]: 400770 : if (indlmn(3,jlmn)/=indlmn(3,ilmn)) cycle
255 : 282510 : jlm = indlmn(4,jlmn)
256 [ + + ]: 282510 : if (ilm<=jlm) then
257 : 177390 : klm_so = jlm*(jlm-1)/2 + ilm
258 : 177390 : sign_so = 1
259 : : else
260 : 105120 : klm_so = ilm*(ilm-1)/2 + jlm
261 : 105120 : sign_so = -1
262 : : end if
263 : 282510 : ls_uu_im = sign_so * ls_ylm_so(2,klm_so,1)
264 : 282510 : ls_ud_re = sign_so * ls_ylm_so(1,klm_so,2)
265 : 282510 : ls_ud_im = sign_so * ls_ylm_so(2,klm_so,2)
266 : : ! up-up: Re(<up|LS|up>)=0, only Im contributes
267 : 282510 : gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1) - ekb_so*ls_uu_im*gx(2,jlmn,ia,1)
268 : 282510 : gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1) + ekb_so*ls_uu_im*gx(1,jlmn,ia,1)
269 : : ! up-dn
270 : : gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1) &
271 : 282510 : & + ekb_so*(ls_ud_re*gx(1,jlmn,ia,2) - ls_ud_im*gx(2,jlmn,ia,2))
272 : : gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1) &
273 : 282510 : & + ekb_so*(ls_ud_re*gx(2,jlmn,ia,2) + ls_ud_im*gx(1,jlmn,ia,2))
274 : : ! dn-up: Re=-ls_ud_re, Im=+ls_ud_im
275 : : gxfac_(1,ilmn,ia,2)=gxfac_(1,ilmn,ia,2) &
276 : 282510 : & + ekb_so*(-ls_ud_re*gx(1,jlmn,ia,1) - ls_ud_im*gx(2,jlmn,ia,1))
277 : : gxfac_(2,ilmn,ia,2)=gxfac_(2,ilmn,ia,2) &
278 : 282510 : & + ekb_so*(-ls_ud_re*gx(2,jlmn,ia,1) + ls_ud_im*gx(1,jlmn,ia,1))
279 : : ! dn-dn: Im=-ls_uu_im
280 : 282510 : gxfac_(1,ilmn,ia,2)=gxfac_(1,ilmn,ia,2) + ekb_so*ls_uu_im*gx(2,jlmn,ia,2)
281 : 1971000 : gxfac_(2,ilmn,ia,2)=gxfac_(2,ilmn,ia,2) - ekb_so*ls_uu_im*gx(1,jlmn,ia,2)
282 : : end do ! jlmn
283 : : end do ! ilmn
284 : : end do ! ia
285 : : end if ! NC+SO
286 : :
287 : : !Accumulate gxfac related to nonlocal operator (PAW)
288 : : !-------------------------------------------------------------------
289 [ + + ]: 38506025 : if (paw_opt==1.or.paw_opt==2.or.paw_opt==4) then
290 : : !Enl is psp strength Dij or (Dij-lambda.Sij)
291 : :
292 : : ! === Diagonal term(s) (up-up, down-down)
293 : :
294 : : ! 1-Enl is real
295 [ + + ]: 29056780 : if (cplex_enl==1) then
296 : : !$OMP PARALLEL &
297 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
298 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,i0lmn,ijlmn,gxi)
299 : : !$OMP DO COLLAPSE(3)
300 [ + + ]: 52622718 : do ispinor=1,nspinor
301 [ + + ]: 87990695 : do ia=1,nincat
302 [ + + ]: 409419499 : do jlmn=1,nlmn
303 : 347740163 : ispinor_index=ispinor+shift
304 : 347740163 : index_enl=atindx1(iatm+ia)
305 : 347740163 : j0lmn=jlmn*(jlmn-1)/2
306 : 347740163 : jjlmn=j0lmn+jlmn
307 : 347740163 : enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
308 [ + + ]: 347740163 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
309 [ + + ]: 1006330769 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
310 [ + + ]: 1006330769 : gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxj(1:cplex)
311 [ + + ]: 2252069564 : do ilmn=1,jlmn-1
312 : 1868961424 : ijlmn=j0lmn+ilmn
313 : 1868961424 : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
314 [ + + ]: 1868961424 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
315 [ + + ]: 5454596461 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
316 [ + + ]: 5454596461 : gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex)
317 : : #if !defined HAVE_OPENMP
318 [ + + ]: 5802336624 : gxfac_(1:cplex,ilmn,ia,ispinor)=gxfac_(1:cplex,ilmn,ia,ispinor)+enl_(1)*gxj(1:cplex)
319 : : #endif
320 : : end do
321 : : #if defined HAVE_OPENMP
322 : : if(jlmn<nlmn) then
323 : : do ilmn=jlmn+1,nlmn
324 : : i0lmn=(ilmn*(ilmn-1)/2)
325 : : ijlmn=i0lmn+jlmn
326 : : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
327 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
328 : : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
329 : : gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex)
330 : : end do
331 : : end if
332 : : #endif
333 : : end do
334 : : end do
335 : : end do
336 : : !$OMP END DO
337 : : !$OMP END PARALLEL
338 : :
339 : : ! 2-Enl is complex ===== D^ss'_ij=D^s's_ji^*
340 : : else
341 [ - + ]: 2745421 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
342 : :
343 [ + + ]: 2745421 : if (nspinortot==1) then ! -------------> NO SPINORS
344 : :
345 : : !$OMP PARALLEL &
346 : : !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,i0lmn,ijlmn,gxi)
347 [ + + ]: 1707526 : do ia=1,nincat
348 : 853763 : index_enl=atindx1(iatm+ia)
349 : : !$OMP DO
350 [ + + ]: 8537630 : do jlmn=1,nlmn
351 : 6830104 : j0lmn=jlmn*(jlmn-1)/2
352 : 6830104 : jjlmn=j0lmn+jlmn
353 : 6830104 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
354 [ - + ]: 6830104 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
355 [ + + ]: 20490312 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,1)
356 : 6830104 : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxj(1)
357 [ + - ]: 6830104 : if (cplex==2) gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxj(2)
358 [ + + ]: 31589231 : do ilmn=1,jlmn-1
359 : 23905364 : ijlmn=j0lmn+ilmn
360 [ + + ]: 71716092 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
361 [ - + ]: 23905364 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
362 [ + + ]: 71716092 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
363 : 23905364 : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxi(1)
364 : 23905364 : gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)-enl_(2)*gxi(1)
365 : : #if !defined HAVE_OPENMP
366 : 23905364 : gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1)+enl_(1)*gxj(1)
367 : 23905364 : gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1)+enl_(2)*gxj(1)
368 : : #endif
369 [ + - ]: 30735468 : if (cplex==2) then
370 : 23905364 : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(2)*gxi(2)
371 : 23905364 : gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxi(2)
372 : : #if !defined HAVE_OPENMP
373 : 23905364 : gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1)-enl_(2)*gxj(2)
374 : 23905364 : gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1)+enl_(1)*gxj(2)
375 : : #endif
376 : : end if
377 : : end do
378 : : #if defined HAVE_OPENMP
379 : : if(jlmn<nlmn) then
380 : : do ilmn=jlmn+1,nlmn
381 : : i0lmn=ilmn*(ilmn-1)/2
382 : : ijlmn=i0lmn+jlmn
383 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
384 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
385 : : gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
386 : : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxi(1)
387 : : gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(2)*gxi(1)
388 : : if (cplex==2) then
389 : : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)-enl_(2)*gxi(2)
390 : : gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxi(2)
391 : : end if
392 : : end do
393 : : end if
394 : : #endif
395 : : end do
396 : : !$OMP END DO
397 : : end do
398 : : !$OMP END PARALLEL
399 : :
400 : : else ! -------------> SPINORIAL CASE
401 : :
402 : : ! === Diagonal term(s) (up-up, down-down)
403 : :
404 : : !$OMP PARALLEL &
405 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
406 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,i0lmn,ijlmn,gxi)
407 [ + + ]: 5614726 : do ispinor=1,nspinor
408 : 3723068 : ispinor_index=ispinor+shift
409 [ + + ]: 10276118 : do ia=1,nincat
410 : 4661392 : index_enl=atindx1(iatm+ia)
411 : : !$OMP DO
412 [ + + ]: 72210884 : do jlmn=1,nlmn
413 : 63826424 : j0lmn=jlmn*(jlmn-1)/2
414 : 63826424 : jjlmn=j0lmn+jlmn
415 : 63826424 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
416 [ + + ]: 63826424 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
417 [ + + ]: 191479272 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
418 : 63826424 : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxj(1)
419 [ + - ]: 63826424 : if (cplex==2) gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxj(2)
420 [ + + ]: 525613256 : do ilmn=1,jlmn-1
421 : 457125440 : ijlmn=j0lmn+ilmn
422 [ + + ]: 1371376320 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
423 [ + + ]: 457125440 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
424 [ + + ]: 1371376320 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
425 : 457125440 : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
426 : 457125440 : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)-enl_(2)*gxi(1)
427 : : #if !defined HAVE_OPENMP
428 : 457125440 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)+enl_(1)*gxj(1)
429 : 457125440 : gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(2)*gxj(1)
430 : : #endif
431 [ + - ]: 520951864 : if (cplex==2) then
432 : 457125440 : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(2)*gxi(2)
433 : 457125440 : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
434 : : #if !defined HAVE_OPENMP
435 : 457125440 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)-enl_(2)*gxj(2)
436 : 457125440 : gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(1)*gxj(2)
437 : : #endif
438 : : end if
439 : : end do
440 : : #if defined HAVE_OPENMP
441 : : if(jlmn<nlmn) then
442 : : do ilmn=jlmn+1,nlmn
443 : : i0lmn=ilmn*(ilmn-1)/2
444 : : ijlmn=i0lmn+jlmn
445 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
446 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
447 : : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
448 : : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
449 : : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(2)*gxi(1)
450 : : if (cplex==2) then
451 : : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)-enl_(2)*gxi(2)
452 : : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
453 : : end if
454 : : end do
455 : : end if
456 : : #endif
457 : : end do
458 : : !$OMP END DO
459 : : end do
460 : : end do
461 : : !$OMP END PARALLEL
462 : : end if !nspinortot
463 : : end if !complex_enl
464 : :
465 : : ! === Off-diagonal term(s) (up-down, down-up)
466 : :
467 : : ! --- No parallelization over spinors ---
468 [ + + + + ]: 29056780 : if (nspinortot==2.and.nspinor==nspinortot) then
469 [ - + ]: 1831410 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
470 [ - + ]: 1831410 : ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex)!")
471 : : !$OMP PARALLEL &
472 : : !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), &
473 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxi,gxj,ilmn,i0lmn,ijlmn)
474 [ + + ]: 5494230 : do ispinor=1,nspinortot
475 : 3662820 : jspinor=3-ispinor
476 [ + + ]: 10052534 : do ia=1,nincat
477 : 4558304 : index_enl=atindx1(iatm+ia)
478 : : !$OMP DO
479 [ + + ]: 70191964 : do jlmn=1,nlmn
480 : 61970840 : j0lmn=jlmn*(jlmn-1)/2
481 : 61970840 : jjlmn=j0lmn+jlmn
482 [ + + ]: 185912520 : enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor )
483 [ + + ]: 185912520 : gxi(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
484 : 61970840 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(1)*gxi(1)
485 : 61970840 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)-enl_(2)*gxi(1)
486 [ + - ]: 61970840 : if (cplex==2) then
487 : 61970840 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(2)*gxi(2)
488 : 61970840 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)+enl_(1)*gxi(2)
489 : : end if
490 : : #if !defined HAVE_OPENMP
491 [ + + ]: 185912520 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,jspinor)
492 : : #endif
493 [ + + ]: 507882120 : do ilmn=1,jlmn-1
494 : 441352976 : ijlmn=j0lmn+ilmn
495 [ + + ]: 1324058928 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
496 [ + + ]: 1324058928 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
497 : 441352976 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(1)*gxi(1)
498 : 441352976 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)-enl_(2)*gxi(1)
499 : : #if !defined HAVE_OPENMP
500 : 441352976 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)+enl_(1)*gxj(1)
501 : 441352976 : gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(2)*gxj(1)
502 : : #endif
503 [ + - ]: 503323816 : if (cplex==2) then
504 : 441352976 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(2)*gxi(2)
505 : 441352976 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)+enl_(1)*gxi(2)
506 : : #if !defined HAVE_OPENMP
507 : 441352976 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)-enl_(2)*gxj(2)
508 : 441352976 : gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(1)*gxj(2)
509 : : #endif
510 : : end if
511 : : end do
512 : : #if defined HAVE_OPENMP
513 : : if(jlmn<nlmn) then
514 : : do ilmn=jlmn+1,nlmn
515 : : i0lmn=ilmn*(ilmn-1)/2
516 : : ijlmn=i0lmn+jlmn
517 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
518 : : gxi(1:cplex)=gx(1:cplex,ilmn,ia,jspinor)
519 : : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
520 : : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(2)*gxi(1)
521 : : if (cplex==2) then
522 : : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)-enl_(2)*gxi(2)
523 : : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
524 : : end if
525 : : end do
526 : : end if
527 : : #endif
528 : : end do
529 : : !$OMP END DO
530 : : end do
531 : : end do
532 : : !$OMP END PARALLEL
533 : :
534 : : ! --- Parallelization over spinors ---
535 [ + + + - ]: 27225370 : else if (nspinortot==2.and.nspinor/=nspinortot) then
536 [ - + ]: 60248 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
537 [ - + ]: 60248 : ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
538 [ - + + - : 361488 : ABI_MALLOC(gxfac_offdiag,(cplex_fac,nlmn,nincat,nspinortot))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
539 : : !$OMP WORKSHARE
540 [ + + + + : 11520424 : gxfac_offdiag(:,:,:,:)=zero
+ + + + ]
541 : : !$OMP END WORKSHARE
542 : 60248 : ispinor_index=mpi_enreg%me_spinor+1
543 : 60248 : jspinor_index=3-ispinor_index
544 [ + + ]: 60248 : if (ispinor_index==1) then
545 : : ijspin=3;jispin=4
546 : : else
547 : 30124 : ijspin=4;jispin=3
548 : : end if
549 : : !$OMP PARALLEL &
550 : : !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,gxi)
551 [ + + ]: 163336 : do ia=1,nincat
552 : 103088 : index_enl=atindx1(iatm+ia)
553 : : !$OMP DO
554 [ + + ]: 2018920 : do jlmn=1,nlmn
555 : 1855584 : j0lmn=jlmn*(jlmn-1)/2
556 [ + + ]: 35359184 : do ilmn=1,nlmn
557 : 33400512 : i0lmn=ilmn*(ilmn-1)/2
558 [ + + ]: 33400512 : if (ilmn<=jlmn) then
559 : 17628048 : ijlmn=j0lmn+ilmn
560 : 17628048 : enl_(1)= enl_ptr(2*ijlmn-1,index_enl,ijspin)
561 : 17628048 : enl_(2)=-enl_ptr(2*ijlmn ,index_enl,ijspin)
562 : : else
563 : 15772464 : jilmn=i0lmn+jlmn
564 : 15772464 : enl_(1)= enl_ptr(2*jilmn-1,index_enl,jispin)
565 : 15772464 : enl_(2)= enl_ptr(2*jilmn ,index_enl,jispin)
566 : : end if
567 [ + + ]: 100201536 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
568 : : gxfac_offdiag(1,jlmn,ia,jspinor_index)= &
569 : 33400512 : & gxfac_offdiag(1,jlmn,ia,jspinor_index)+enl_(1)*gxi(1)
570 : : gxfac_offdiag(2,jlmn,ia,jspinor_index)= &
571 : 33400512 : & gxfac_offdiag(2,jlmn,ia,jspinor_index)+enl_(2)*gxi(1)
572 [ + - ]: 35256096 : if (cplex==2) then
573 : : gxfac_offdiag(1,jlmn,ia,jspinor_index)= &
574 : 33400512 : & gxfac_offdiag(1,jlmn,ia,jspinor_index)-enl_(2)*gxi(2)
575 : : gxfac_offdiag(2,jlmn,ia,jspinor_index)= &
576 : 33400512 : & gxfac_offdiag(2,jlmn,ia,jspinor_index)+enl_(1)*gxi(2)
577 : : end if
578 : : end do !ilmn
579 : : end do !jlmn
580 : : !$OMP END DO
581 : : end do !iat
582 : : !$OMP END PARALLEL
583 : 60248 : call xmpi_sum(gxfac_offdiag,mpi_enreg%comm_spinor,ierr)
584 [ + + + + : 5730088 : gxfac_(:,:,:,1)=gxfac_(:,:,:,1)+gxfac_offdiag(:,:,:,ispinor_index)
+ + ]
585 [ - + ]: 120496 : ABI_FREE(gxfac_offdiag)
586 : : end if
587 : :
588 : : end if !paw_opt
589 : :
590 : : !Accumulate dgxdtfac related to nonlocal operator (Norm-conserving)
591 : : !-------------------------------------------------------------------
592 [ + + + + ]: 38506025 : if (optder>=1.and.paw_opt==0) then
593 : : !Enl is E(Kleinman-Bylander)
594 [ - + ]: 2665549 : ABI_CHECK(cplex_enl==1,"BUG: invalid cplex_enl/=1!")
595 [ - + ]: 2665549 : ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!")
596 : : !$OMP PARALLEL &
597 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_,mu)
598 [ + + ]: 5332842 : do ispinor=1,nspinor
599 : 2667293 : ispinor_index = ispinor + shift
600 [ + + ]: 8909133 : do ia=1,nincat
601 : : !$OMP DO
602 [ + + ]: 39681661 : do ilmn=1,nlmn
603 [ + + ]: 33438077 : if (indlmn(6,ilmn)==2) cycle ! NC+SO: SO projectors handled separately below
604 : 33418893 : iln=indlmn(5,ilmn)
605 : 33418893 : enl_(1)=enl_ptr(iln,itypat,ispinor_index)
606 [ + + ]: 78514750 : do mu=1,ndgxdtfac
607 [ + + ]: 157996775 : dgxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=enl_(1)*dgxdt(1:cplex,mu,ilmn,ia,ispinor)
608 : : end do
609 : : end do
610 : : !$OMP END DO
611 : : end do
612 : : end do
613 : : !$OMP END PARALLEL
614 : : end if
615 : :
616 : : !--- NC+SO: real-Ylm L.S coupling for first derivatives ---
617 [ + + + + ]: 38506025 : if (optder>=1.and.lmax_so > 0) then
618 [ + + ]: 1744 : do ia=1,nincat
619 [ + + ]: 23544 : do ilmn=1,nlmn
620 [ + + ]: 21800 : if (indlmn(6,ilmn)/=2) cycle
621 : 9592 : iln = indlmn(5,ilmn)
622 : 9592 : ekb_so = enl_ptr(iln,itypat,1)
623 [ - + ]: 9592 : if (abs(ekb_so)<tol16) cycle
624 : 9592 : ll_so = indlmn(1,ilmn)
625 : 9592 : ilm = indlmn(4,ilmn)
626 [ + + ]: 250264 : do jlmn=1,nlmn
627 [ + + ]: 239800 : if (indlmn(6,jlmn)/=2) cycle
628 [ + + ]: 105512 : if (indlmn(1,jlmn)/=ll_so) cycle
629 [ + + ]: 53192 : if (indlmn(3,jlmn)/=indlmn(3,ilmn)) cycle
630 : 37496 : jlm = indlmn(4,jlmn)
631 [ + + ]: 37496 : if (ilm<=jlm) then
632 : 23544 : klm_so = jlm*(jlm-1)/2 + ilm
633 : 23544 : sign_so = 1
634 : : else
635 : 13952 : klm_so = ilm*(ilm-1)/2 + jlm
636 : 13952 : sign_so = -1
637 : : end if
638 : 37496 : ls_uu_im = sign_so * ls_ylm_so(2,klm_so,1)
639 : 37496 : ls_ud_re = sign_so * ls_ylm_so(1,klm_so,2)
640 : 37496 : ls_ud_im = sign_so * ls_ylm_so(2,klm_so,2)
641 [ + + ]: 99544 : do mu=1,ndgxdtfac
642 : 40248 : dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1) - ekb_so*ls_uu_im*dgxdt(2,mu,jlmn,ia,1)
643 : 40248 : dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1) + ekb_so*ls_uu_im*dgxdt(1,mu,jlmn,ia,1)
644 : : dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1) &
645 : 40248 : & + ekb_so*(ls_ud_re*dgxdt(1,mu,jlmn,ia,2) - ls_ud_im*dgxdt(2,mu,jlmn,ia,2))
646 : : dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1) &
647 : 40248 : & + ekb_so*(ls_ud_re*dgxdt(2,mu,jlmn,ia,2) + ls_ud_im*dgxdt(1,mu,jlmn,ia,2))
648 : : dgxdtfac_(1,mu,ilmn,ia,2)=dgxdtfac_(1,mu,ilmn,ia,2) &
649 : 40248 : & + ekb_so*(-ls_ud_re*dgxdt(1,mu,jlmn,ia,1) - ls_ud_im*dgxdt(2,mu,jlmn,ia,1))
650 : : dgxdtfac_(2,mu,ilmn,ia,2)=dgxdtfac_(2,mu,ilmn,ia,2) &
651 : 40248 : & + ekb_so*(-ls_ud_re*dgxdt(2,mu,jlmn,ia,1) + ls_ud_im*dgxdt(1,mu,jlmn,ia,1))
652 : 40248 : dgxdtfac_(1,mu,ilmn,ia,2)=dgxdtfac_(1,mu,ilmn,ia,2) + ekb_so*ls_uu_im*dgxdt(2,mu,jlmn,ia,2)
653 : 280048 : dgxdtfac_(2,mu,ilmn,ia,2)=dgxdtfac_(2,mu,ilmn,ia,2) - ekb_so*ls_uu_im*dgxdt(1,mu,jlmn,ia,2)
654 : : end do ! mu
655 : : end do ! jlmn
656 : : end do ! ilmn
657 : : end do ! ia
658 : : end if ! NC+SO dgxdtfac_
659 : :
660 : : !Accumulate dgxdtfac related to nonlocal operator (PAW)
661 : : !-------------------------------------------------------------------
662 [ + + + + ]: 38506025 : if (optder>=1.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then
663 : : !Enl is psp strength Dij or (Dij-lambda.Sij)
664 : :
665 : : ! === Diagonal term(s) (up-up, down-down)
666 : :
667 : : ! 1-Enl is real
668 [ + + ]: 1458897 : if (cplex_enl==1) then
669 : : !$OMP PARALLEL &
670 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
671 [ - + + - : 5063196 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
- + + - +
- + - - +
- + ]
672 [ + + ]: 2531598 : do ispinor=1,nspinor
673 : 1265799 : ispinor_index=ispinor+shift
674 [ + + ]: 3995348 : do ia=1,nincat
675 : 1463750 : index_enl=atindx1(iatm+ia)
676 : : !$OMP DO
677 [ + + ]: 15593735 : do jlmn=1,nlmn
678 : 12864186 : j0lmn=jlmn*(jlmn-1)/2
679 : 12864186 : jjlmn=j0lmn+jlmn
680 : 12864186 : enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
681 [ + + ]: 12864186 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
682 [ + + ]: 26397914 : do mu=1,ndgxdtfac
683 [ + + ]: 40601184 : gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
684 [ + + ]: 53465370 : dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
685 : : end do
686 [ + + ]: 70489832 : do ilmn=1,jlmn-1
687 : 56161896 : ijlmn=j0lmn+ilmn
688 : 56161896 : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
689 [ + + ]: 56161896 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
690 [ + + ]: 128261486 : do mu=1,ndgxdtfac
691 [ + + ]: 177706212 : gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
692 [ + + ]: 177706212 : dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
693 : : #if !defined HAVE_OPENMP
694 [ + + ]: 233868108 : dgxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
695 : : #endif
696 : : end do
697 : : end do
698 : : #if defined HAVE_OPENMP
699 : : if(jlmn<nlmn) then
700 : : do ilmn=jlmn+1,nlmn
701 : : i0lmn=ilmn*(ilmn-1)/2
702 : : ijlmn=i0lmn+jlmn
703 : : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
704 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
705 : : do mu=1,ndgxdtfac
706 : : gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
707 : : dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
708 : : end do
709 : : end do
710 : : end if
711 : : #endif
712 : : end do
713 : : !$OMP END DO
714 : : end do
715 : : end do
716 [ - + ]: 1265799 : ABI_FREE(gxfj)
717 : : !$OMP END PARALLEL
718 : :
719 : : ! 2-Enl is complex ===== D^ss'_ij=D^s's_ji^*
720 : : else
721 [ - + ]: 193098 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
722 : :
723 [ + + ]: 193098 : if (nspinortot==1) then ! -------------> NO SPINORS
724 : :
725 : : !$OMP PARALLEL &
726 : : !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
727 [ - + + - : 490224 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
- + + - +
- + - - +
- + ]
728 [ + + ]: 245112 : do ia=1,nincat
729 : 122556 : index_enl=atindx1(iatm+ia)
730 : : !$OMP DO
731 [ + + ]: 1225560 : do jlmn=1,nlmn
732 : 980448 : j0lmn=jlmn*(jlmn-1)/2
733 : 980448 : jjlmn=j0lmn+jlmn
734 : 980448 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
735 [ - + ]: 980448 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
736 [ + + ]: 1960896 : do mu=1,ndgxdtfac
737 [ - + ]: 980448 : if(cplex_dgxdt(mu)==2)then
738 : 0 : cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,1)
739 : : else
740 [ + + ]: 2941344 : cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,1)
741 : : end if
742 : 980448 : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu)
743 [ + - ]: 1960896 : if (cplex_==2) dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu)
744 : : end do
745 [ + + ]: 4534572 : do ilmn=1,jlmn-1
746 : 3431568 : ijlmn=j0lmn+ilmn
747 [ + + ]: 10294704 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
748 [ - + ]: 3431568 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
749 [ + + ]: 7843584 : do mu=1,ndgxdtfac
750 [ - + ]: 3431568 : if(cplex_dgxdt(mu)==2)then
751 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
752 : : else
753 [ + + ]: 10294704 : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
754 : : end if
755 : 3431568 : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
756 : 3431568 : dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)-enl_(2)*gxfi(1)
757 : : #if !defined HAVE_OPENMP
758 : 3431568 : dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu)
759 : 3431568 : dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu)
760 : : #endif
761 [ + - ]: 6863136 : if (cplex_==2) then
762 : 3431568 : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(2)*gxfi(2)
763 : 3431568 : dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
764 : : #if !defined HAVE_OPENMP
765 : 3431568 : dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu)
766 : 3431568 : dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1)+enl_(1)*gxfj(2,mu)
767 : : #endif
768 : : end if
769 : : end do
770 : : end do
771 : : #if defined HAVE_OPENMP
772 : : if(jlmn<nlmn) then
773 : : do ilmn=jlmn+1,nlmn
774 : : i0lmn=ilmn*(ilmn-1)/2
775 : : ijlmn=i0lmn+jlmn
776 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
777 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
778 : : do mu=1,ndgxdtfac
779 : : if(cplex_dgxdt(mu)==2)then
780 : : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
781 : : else
782 : : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
783 : : end if
784 : : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
785 : : dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1)
786 : : if (cplex_==2) then
787 : : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2)
788 : : dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
789 : : end if
790 : : end do
791 : : end do
792 : : end if
793 : : #endif
794 : : end do
795 : : !$OMP END DO
796 : : end do
797 [ - + ]: 122556 : ABI_FREE(gxfj)
798 : : !$OMP END PARALLEL
799 : :
800 : : else ! -------------> SPINORIAL CASE
801 : :
802 : : ! === Diagonal term(s) (up-up, down-down)
803 : :
804 : : !$OMP PARALLEL &
805 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
806 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
807 [ - + + - : 282168 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
- + + - +
- + - - +
- + ]
808 [ + + ]: 211626 : do ispinor=1,nspinor
809 : 141084 : ispinor_index = ispinor + shift
810 [ + + ]: 361698 : do ia=1,nincat
811 : 150072 : index_enl=atindx1(iatm+ia)
812 : : !$OMP DO
813 [ + + ]: 1581612 : do jlmn=1,nlmn
814 : 1290456 : j0lmn=jlmn*(jlmn-1)/2
815 : 1290456 : jjlmn=j0lmn+jlmn
816 : 1290456 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
817 [ - + ]: 1290456 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
818 [ + + ]: 2594952 : do mu=1,ndgxdtfac
819 [ - + ]: 1304496 : if(cplex_dgxdt(mu)==2)then
820 : 0 : cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,ispinor)
821 : : else
822 [ + + ]: 3913488 : cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
823 : : end if
824 : 1304496 : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
825 [ + - ]: 2594952 : if (cplex_==2) dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
826 : : end do
827 [ + + ]: 6541344 : do ilmn=1,jlmn-1
828 : 5100816 : ijlmn=j0lmn+ilmn
829 [ + + ]: 15302448 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
830 [ - + ]: 5100816 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
831 [ + + ]: 11576328 : do mu=1,ndgxdtfac
832 [ - + ]: 5185056 : if(cplex_dgxdt(mu)==2)then
833 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor)
834 : : else
835 [ + + ]: 15555168 : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
836 : : end if
837 : 5185056 : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
838 : 5185056 : dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(1)
839 : : #if !defined HAVE_OPENMP
840 : 5185056 : dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
841 : 5185056 : dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
842 : : #endif
843 [ + - ]: 10285872 : if (cplex_==2) then
844 : 5185056 : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(2)
845 : 5185056 : dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
846 : : #if !defined HAVE_OPENMP
847 : 5185056 : dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
848 : 5185056 : dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
849 : : #endif
850 : : end if
851 : : end do
852 : : end do
853 : : #if defined HAVE_OPENMP
854 : : if(jlmn<nlmn) then
855 : : do ilmn=jlmn+1,nlmn
856 : : i0lmn=ilmn*(ilmn-1)/2
857 : : ijlmn=i0lmn+jlmn
858 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
859 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
860 : : do mu=1,ndgxdtfac
861 : : if(cplex_dgxdt(mu)==2)then
862 : : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor)
863 : : else
864 : : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
865 : : end if
866 : : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
867 : : dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
868 : : if (cplex_==2) then
869 : : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
870 : : dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
871 : : end if
872 : : end do
873 : : end do
874 : : end if
875 : : #endif
876 : : end do
877 : : !$OMP END DO
878 : : end do
879 : : end do
880 [ - + ]: 70542 : ABI_FREE(gxfj)
881 : : !$OMP END PARALLEL
882 : : end if !nspinortot
883 : : end if !complex
884 : :
885 : : ! === Off-diagonal term(s) (up-down, down-up)
886 : :
887 : : ! --- No parallelization over spinors ---
888 [ + + + - ]: 1458897 : if (nspinortot==2.and.nspinor==nspinortot) then
889 [ - + ]: 70542 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
890 [ - + ]: 70542 : ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
891 : : !$OMP PARALLEL &
892 : : !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), &
893 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfi,gxfj,ilmn,i0lmn,ijlmn)
894 [ - + + - : 282168 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
- + + - +
- + - - +
- + ]
895 [ + + ]: 211626 : do ispinor=1,nspinor
896 : 141084 : jspinor=3-ispinor
897 [ + + ]: 361698 : do ia=1,nincat
898 : 150072 : index_enl=atindx1(iatm+ia)
899 : : !$OMP DO
900 [ + + ]: 1581612 : do jlmn=1,nlmn
901 : 1290456 : j0lmn=jlmn*(jlmn-1)/2
902 : 1290456 : jjlmn=j0lmn+jlmn
903 [ + + ]: 3871368 : enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor)
904 [ + + ]: 2594952 : do mu=1,ndgxdtfac
905 [ - + ]: 1304496 : if(cplex_dgxdt(mu)==2)then
906 : 0 : cplex_ = 2 ;
907 : 0 : gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,jlmn,ia,ispinor)
908 : 0 : gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,jspinor)
909 : : else
910 : 1304496 : cplex_ = cplex ;
911 [ + + ]: 3913488 : gxfi(1:cplex) =dgxdt(1:cplex,mu,jlmn,ia,ispinor)
912 [ + + ]: 3913488 : gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,jspinor)
913 : : end if
914 : 1304496 : dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
915 : 1304496 : dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
916 [ + - ]: 2594952 : if (cplex_==2) then
917 : 1304496 : dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
918 : 1304496 : dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
919 : : end if
920 : : end do
921 [ + + ]: 6541344 : do ilmn=1,jlmn-1
922 : 5100816 : ijlmn=j0lmn+ilmn
923 [ + + ]: 15302448 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
924 [ + + ]: 11576328 : do mu=1,ndgxdtfac
925 [ - + ]: 5185056 : if(cplex_dgxdt(mu)==2)then
926 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor)
927 : : else
928 [ + + ]: 15555168 : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
929 : : end if
930 : 5185056 : dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
931 : 5185056 : dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
932 : : #if !defined HAVE_OPENMP
933 : 5185056 : dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
934 : 5185056 : dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
935 : : #endif
936 [ + - ]: 10285872 : if (cplex_==2) then
937 : 5185056 : dgxdtfac_(1,mu,jlmn,ia,jspinor)=dgxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
938 : 5185056 : dgxdtfac_(2,mu,jlmn,ia,jspinor)=dgxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
939 : : #if !defined HAVE_OPENMP
940 : 5185056 : dgxdtfac_(1,mu,ilmn,ia,ispinor)=dgxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
941 : 5185056 : dgxdtfac_(2,mu,ilmn,ia,ispinor)=dgxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
942 : : #endif
943 : : end if
944 : : end do !mu
945 : : end do !ilmn
946 : : #if defined HAVE_OPENMP
947 : : if(jlmn<nlmn) then
948 : : do ilmn=jlmn+1,nlmn
949 : : i0lmn=ilmn*(ilmn-1)/2
950 : : ijlmn=i0lmn+jlmn
951 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
952 : : do mu=1,ndgxdtfac
953 : : if(cplex_dgxdt(mu)==2)then
954 : : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,jspinor)
955 : : else
956 : : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,jspinor)
957 : : end if
958 : : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
959 : : dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
960 : : if (cplex_==2) then
961 : : dgxdtfac_(1,mu,jlmn,ia,ispinor)=dgxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
962 : : dgxdtfac_(2,mu,jlmn,ia,ispinor)=dgxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
963 : : end if
964 : : end do !mu
965 : : end do !ilmn
966 : : end if
967 : : #endif
968 : : end do !jmln
969 : : !$OMP END DO
970 : : end do !ia
971 : : end do !ispinor
972 [ - + ]: 70542 : ABI_FREE(gxfj)
973 : : !$OMP END PARALLEL
974 : :
975 : : ! --- Parallelization over spinors ---
976 [ - + - - ]: 1388355 : else if (nspinortot==2.and.nspinor/=nspinortot) then
977 [ # # ]: 0 : ABI_CHECK(cplex_enl==2,"BUG: opernlc_ylm: invalid cplex_enl/=2!")
978 [ # # ]: 0 : ABI_CHECK(cplex_fac==2,"BUG: opernlc_ylm: invalid cplex_fac/=2!")
979 [ # # # # : 0 : ABI_MALLOC(dgxdtfac_offdiag,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinortot))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
980 : : !$OMP PARALLEL &
981 : : !$OMP PRIVATE(ia,index_enl), &
982 : : !$OMP PRIVATE(jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,mu,gxfi)
983 : : !$OMP WORKSHARE
984 [ # # # # : 0 : dgxdtfac_offdiag(:,:,:,:,:)=zero
# # # # #
# ]
985 : : !$OMP END WORKSHARE
986 : : !$OMP SINGLE
987 : 0 : ispinor_index=mpi_enreg%me_spinor+1
988 : 0 : jspinor_index=3-ispinor_index
989 [ # # ]: 0 : if (ispinor_index==1) then
990 : : ijspin=3;jispin=4
991 : : else
992 : 0 : ijspin=4;jispin=3
993 : : end if
994 : : !$OMP END SINGLE
995 [ # # ]: 0 : do ia=1,nincat
996 : 0 : index_enl=atindx1(iatm+ia)
997 : : !$OMP DO
998 [ # # ]: 0 : do jlmn=1,nlmn
999 : 0 : j0lmn=jlmn*(jlmn-1)/2
1000 [ # # ]: 0 : do ilmn=1,nlmn
1001 : 0 : i0lmn=ilmn*(ilmn-1)/2
1002 [ # # ]: 0 : if (ilmn<=jlmn) then
1003 : 0 : ijlmn=j0lmn+ilmn
1004 : 0 : enl_(1)= enl_ptr(2*ijlmn-1,index_enl,ijspin)
1005 : 0 : enl_(2)=-enl_ptr(2*ijlmn ,index_enl,ijspin)
1006 : : else
1007 : 0 : jilmn=i0lmn+jlmn
1008 : 0 : enl_(1)= enl_ptr(2*jilmn-1,index_enl,jispin)
1009 : 0 : enl_(2)= enl_ptr(2*jilmn ,index_enl,jispin)
1010 : : end if
1011 [ # # ]: 0 : do mu=1,ndgxdtfac
1012 [ # # ]: 0 : if(cplex_dgxdt(mu)==2)then
1013 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
1014 : : else
1015 [ # # ]: 0 : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
1016 : : end if
1017 : : dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
1018 : 0 : & dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(1)
1019 : : dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
1020 : 0 : & dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(2)*gxfi(1)
1021 [ # # ]: 0 : if (cplex_==2) then
1022 : : dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
1023 : 0 : & dgxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)-enl_(2)*gxfi(2)
1024 : : dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
1025 : 0 : & dgxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(2)
1026 : : end if
1027 : : end do
1028 : : end do !ilmn
1029 : : end do !jlmn
1030 : : !$OMP END DO
1031 : : end do !iat
1032 : : !$OMP SINGLE
1033 : 0 : call xmpi_sum(dgxdtfac_offdiag,mpi_enreg%comm_spinor,ierr)
1034 : : !$OMP END SINGLE
1035 : : !$OMP WORKSHARE
1036 [ # # # # : 0 : dgxdtfac_(:,:,:,:,1)=dgxdtfac_(:,:,:,:,1)+dgxdtfac_offdiag(:,:,:,:,ispinor_index)
# # # # ]
1037 : : !$OMP END WORKSHARE
1038 : : !$OMP END PARALLEL
1039 [ # # ]: 0 : ABI_FREE(dgxdtfac_offdiag)
1040 : : end if !nspinortot
1041 : :
1042 : : end if ! pawopt & optder
1043 : :
1044 : : !Accumulate d2gxdtfac related to nonlocal operator (Norm-conserving)
1045 : : !-------------------------------------------------------------------
1046 [ + + + + ]: 38506025 : if (optder==2.and.paw_opt==0) then
1047 : : !Enl is E(Kleinman-Bylander)
1048 : : !$OMP PARALLEL &
1049 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_,mu)
1050 [ + + ]: 898032 : do ispinor=1,nspinor
1051 : 449016 : ispinor_index = ispinor + shift
1052 [ + + ]: 1742067 : do ia=1,nincat
1053 : : !$OMP DO
1054 [ + + ]: 8682588 : do ilmn=1,nlmn
1055 [ + - ]: 7389537 : if (indlmn(6,ilmn)==2) cycle ! NC+SO: SO projectors handled separately below
1056 : 7389537 : iln=indlmn(5,ilmn)
1057 : 7389537 : enl_(1)=enl_ptr(iln,itypat,ispinor_index)
1058 [ + + ]: 28562373 : do mu=1,nd2gxdtfac
1059 [ + + ]: 68375940 : d2gxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=enl_(1)*d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1060 : : end do
1061 : : end do
1062 : : !$OMP END DO
1063 : : end do
1064 : : end do
1065 : : !$OMP END PARALLEL
1066 : : end if
1067 : :
1068 : : ! NC+SO: real-Ylm L.S coupling for second derivatives (elastic tensor) ---
1069 [ + + - + ]: 38506025 : if (optder==2.and.lmax_so > 0) then
1070 [ # # ]: 0 : do ia=1,nincat
1071 [ # # ]: 0 : do ilmn=1,nlmn
1072 [ # # ]: 0 : if (indlmn(6,ilmn)/=2) cycle
1073 : 0 : iln = indlmn(5,ilmn)
1074 : 0 : ekb_so = enl_ptr(iln,itypat,1)
1075 [ # # ]: 0 : if (abs(ekb_so)<tol16) cycle
1076 : 0 : ll_so = indlmn(1,ilmn)
1077 : 0 : ilm = indlmn(4,ilmn)
1078 [ # # ]: 0 : do jlmn=1,nlmn
1079 [ # # ]: 0 : if (indlmn(6,jlmn)/=2) cycle
1080 [ # # ]: 0 : if (indlmn(1,jlmn)/=ll_so) cycle
1081 [ # # ]: 0 : if (indlmn(3,jlmn)/=indlmn(3,ilmn)) cycle
1082 : 0 : jlm = indlmn(4,jlmn)
1083 [ # # ]: 0 : if (ilm<=jlm) then
1084 : 0 : klm_so = jlm*(jlm-1)/2 + ilm
1085 : 0 : sign_so = 1
1086 : : else
1087 : 0 : klm_so = ilm*(ilm-1)/2 + jlm
1088 : 0 : sign_so = -1
1089 : : end if
1090 : 0 : ls_uu_im = sign_so * ls_ylm_so(2,klm_so,1)
1091 : 0 : ls_ud_re = sign_so * ls_ylm_so(1,klm_so,2)
1092 : 0 : ls_ud_im = sign_so * ls_ylm_so(2,klm_so,2)
1093 [ # # ]: 0 : do mu=1,nd2gxdtfac
1094 : : ! up-up
1095 : 0 : d2gxdtfac_(1,mu,ilmn,ia,1)=d2gxdtfac_(1,mu,ilmn,ia,1) - ekb_so*ls_uu_im*d2gxdt(2,mu,jlmn,ia,1)
1096 : 0 : d2gxdtfac_(2,mu,ilmn,ia,1)=d2gxdtfac_(2,mu,ilmn,ia,1) + ekb_so*ls_uu_im*d2gxdt(1,mu,jlmn,ia,1)
1097 : : ! up-dn
1098 : : d2gxdtfac_(1,mu,ilmn,ia,1)=d2gxdtfac_(1,mu,ilmn,ia,1) &
1099 : 0 : & + ekb_so*(ls_ud_re*d2gxdt(1,mu,jlmn,ia,2) - ls_ud_im*d2gxdt(2,mu,jlmn,ia,2))
1100 : : d2gxdtfac_(2,mu,ilmn,ia,1)=d2gxdtfac_(2,mu,ilmn,ia,1) &
1101 : 0 : & + ekb_so*(ls_ud_re*d2gxdt(2,mu,jlmn,ia,2) + ls_ud_im*d2gxdt(1,mu,jlmn,ia,2))
1102 : : ! dn-up
1103 : : d2gxdtfac_(1,mu,ilmn,ia,2)=d2gxdtfac_(1,mu,ilmn,ia,2) &
1104 : 0 : & + ekb_so*(-ls_ud_re*d2gxdt(1,mu,jlmn,ia,1) - ls_ud_im*d2gxdt(2,mu,jlmn,ia,1))
1105 : : d2gxdtfac_(2,mu,ilmn,ia,2)=d2gxdtfac_(2,mu,ilmn,ia,2) &
1106 : 0 : & + ekb_so*(-ls_ud_re*d2gxdt(2,mu,jlmn,ia,1) + ls_ud_im*d2gxdt(1,mu,jlmn,ia,1))
1107 : : ! dn-dn
1108 : 0 : d2gxdtfac_(1,mu,ilmn,ia,2)=d2gxdtfac_(1,mu,ilmn,ia,2) + ekb_so*ls_uu_im*d2gxdt(2,mu,jlmn,ia,2)
1109 : 0 : d2gxdtfac_(2,mu,ilmn,ia,2)=d2gxdtfac_(2,mu,ilmn,ia,2) - ekb_so*ls_uu_im*d2gxdt(1,mu,jlmn,ia,2)
1110 : : end do ! mu
1111 : : end do ! jlmn
1112 : : end do ! ilmn
1113 : : end do ! ia
1114 : : end if ! NC+SO d2gxdtfac
1115 : :
1116 [ + + ]: 38506025 : if (lmax_so > 0) then
1117 [ - + ]: 6570 : ABI_FREE(ls_ylm_so)
1118 : : end if
1119 : :
1120 : : DBG_EXIT("COLL")
1121 : :
1122 : : !Accumulate d2gxdtfac related to nonlocal operator (PAW)
1123 : : !-------------------------------------------------------------------
1124 [ + + + + ]: 38506025 : if (optder==2.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then
1125 : : !Enl is psp strength Dij or (Dij-lambda.Sij)
1126 : :
1127 : : ! === Diagonal term(s) (up-up, down-down)
1128 : :
1129 : : ! 1-Enl is real
1130 [ + + ]: 4173 : if (cplex_enl==1) then
1131 : : !$OMP PARALLEL &
1132 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
1133 [ - + + - : 15612 : ABI_MALLOC(gxfj,(cplex,nd2gxdtfac))
- + + - +
- + - - +
- + ]
1134 [ + + ]: 7806 : do ispinor=1,nspinor
1135 : 3903 : ispinor_index=ispinor+shift
1136 [ + + ]: 11772 : do ia=1,nincat
1137 : 3966 : index_enl=atindx1(iatm+ia)
1138 : : !$OMP DO
1139 [ + + ]: 40227 : do jlmn=1,nlmn
1140 : 32358 : j0lmn=jlmn*(jlmn-1)/2
1141 : 32358 : jjlmn=j0lmn+jlmn
1142 : 32358 : enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
1143 [ - + ]: 32358 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
1144 [ + + ]: 64716 : do mu=1,nd2gxdtfac
1145 [ + + ]: 97074 : gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1146 [ + + ]: 129432 : d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
1147 : : end do
1148 [ + + ]: 153672 : do ilmn=1,jlmn-1
1149 : 117348 : ijlmn=j0lmn+ilmn
1150 : 117348 : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
1151 [ - + ]: 117348 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1152 [ + + ]: 267054 : do mu=1,nd2gxdtfac
1153 [ + + ]: 352044 : gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1154 [ + + ]: 352044 : d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
1155 : : #if !defined HAVE_OPENMP
1156 [ + + ]: 469392 : d2gxdtfac_(1:cplex,mu,ilmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
1157 : : #endif
1158 : : end do
1159 : : end do
1160 : : #if defined HAVE_OPENMP
1161 : : if(jlmn<nlmn) then
1162 : : do ilmn=jlmn+1,nlmn
1163 : : i0lmn=ilmn*(ilmn-1)/2
1164 : : ijlmn=i0lmn+jlmn
1165 : : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
1166 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1167 : : do mu=1,nd2gxdtfac
1168 : : gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1169 : : d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex)
1170 : : end do
1171 : : end do
1172 : : end if
1173 : : #endif
1174 : : end do
1175 : : !$OMP END DO
1176 : : end do
1177 : : end do
1178 [ - + ]: 3903 : ABI_FREE(gxfj)
1179 : : !$OMP END PARALLEL
1180 : :
1181 : : ! 2-Enl is complex ===== D^ss'_ij=D^s's_ji^*
1182 : : else
1183 [ - + ]: 270 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
1184 : :
1185 [ - + ]: 270 : if (nspinortot==1) then ! -------------> NO SPINORS
1186 : :
1187 : : !$OMP PARALLEL &
1188 : : !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
1189 [ # # # # : 0 : ABI_MALLOC(gxfj,(cplex,nd2gxdtfac))
# # # # #
# # # # #
# # ]
1190 [ # # ]: 0 : do ia=1,nincat
1191 : 0 : index_enl=atindx1(iatm+ia)
1192 : : !$OMP DO
1193 [ # # ]: 0 : do jlmn=1,nlmn
1194 : 0 : j0lmn=jlmn*(jlmn-1)/2
1195 : 0 : jjlmn=j0lmn+jlmn
1196 : 0 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
1197 [ # # ]: 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
1198 [ # # ]: 0 : do mu=1,nd2gxdtfac
1199 [ # # ]: 0 : if(cplex_d2gxdt(mu)==2)then
1200 : 0 : cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,1)
1201 : : else
1202 [ # # ]: 0 : cplex_ = cplex ; gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,1)
1203 : : end if
1204 : 0 : d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu)
1205 [ # # ]: 0 : if (cplex_==2) d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu)
1206 : : end do
1207 [ # # ]: 0 : do ilmn=1,jlmn-1
1208 : 0 : ijlmn=j0lmn+ilmn
1209 [ # # ]: 0 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
1210 [ # # ]: 0 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1211 [ # # ]: 0 : do mu=1,nd2gxdtfac
1212 [ # # ]: 0 : if(cplex_d2gxdt(mu)==2)then
1213 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1)
1214 : : else
1215 [ # # ]: 0 : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1)
1216 : : end if
1217 : 0 : d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
1218 : 0 : d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)-enl_(2)*gxfi(1)
1219 : : #if !defined HAVE_OPENMP
1220 : 0 : d2gxdtfac_(1,mu,ilmn,ia,1)=d2gxdtfac_(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu)
1221 : 0 : d2gxdtfac_(2,mu,ilmn,ia,1)=d2gxdtfac_(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu)
1222 : : #endif
1223 [ # # ]: 0 : if (cplex_==2) then
1224 : 0 : d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(2)*gxfi(2)
1225 : 0 : d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
1226 : : #if !defined HAVE_OPENMP
1227 : 0 : d2gxdtfac_(1,mu,ilmn,ia,1)=d2gxdtfac_(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu)
1228 : 0 : d2gxdtfac_(2,mu,ilmn,ia,1)=d2gxdtfac_(2,mu,ilmn,ia,1)+enl_(1)*gxfj(2,mu)
1229 : : #endif
1230 : : end if
1231 : : end do
1232 : : end do
1233 : : #if defined HAVE_OPENMP
1234 : : if(jlmn<nlmn) then
1235 : : do ilmn=jlmn+1,nlmn
1236 : : i0lmn=ilmn*(ilmn-1)/2
1237 : : ijlmn=i0lmn+jlmn
1238 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
1239 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1240 : : do mu=1,nd2gxdtfac
1241 : : if(cplex_d2gxdt(mu)==2)then
1242 : : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1)
1243 : : else
1244 : : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1)
1245 : : end if
1246 : : d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
1247 : : d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1)
1248 : : if (cplex_==2) then
1249 : : d2gxdtfac_(1,mu,jlmn,ia,1)=d2gxdtfac_(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2)
1250 : : d2gxdtfac_(2,mu,jlmn,ia,1)=d2gxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
1251 : : end if
1252 : : end do
1253 : : end do
1254 : : end if
1255 : : #endif
1256 : : end do
1257 : : !$OMP END DO
1258 : : end do
1259 [ # # ]: 0 : ABI_FREE(gxfj)
1260 : : !$OMP END PARALLEL
1261 : :
1262 : : else ! -------------> SPINORIAL CASE
1263 : :
1264 : : ! === Diagonal term(s) (up-up, down-down)
1265 : :
1266 : : !$OMP PARALLEL &
1267 : : !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), &
1268 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
1269 [ - + + - : 1080 : ABI_MALLOC(gxfj,(cplex,nd2gxdtfac))
- + + - +
- + - - +
- + ]
1270 [ + + ]: 810 : do ispinor=1,nspinor
1271 : 540 : ispinor_index = ispinor + shift
1272 [ + + ]: 1890 : do ia=1,nincat
1273 : 1080 : index_enl=atindx1(iatm+ia)
1274 : : !$OMP DO
1275 [ + + ]: 15660 : do jlmn=1,nlmn
1276 : 14040 : j0lmn=jlmn*(jlmn-1)/2
1277 : 14040 : jjlmn=j0lmn+jlmn
1278 : 14040 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
1279 [ - + ]: 14040 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
1280 [ + + ]: 28080 : do mu=1,nd2gxdtfac
1281 [ - + ]: 14040 : if(cplex_d2gxdt(mu)==2)then
1282 : 0 : cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,ispinor)
1283 : : else
1284 [ + + ]: 42120 : cplex_ = cplex ; gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1285 : : end if
1286 : 14040 : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
1287 [ + - ]: 28080 : if (cplex_==2) d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
1288 : : end do
1289 [ + + ]: 99360 : do ilmn=1,jlmn-1
1290 : 84240 : ijlmn=j0lmn+ilmn
1291 [ + + ]: 252720 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
1292 [ - + ]: 84240 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1293 [ + + ]: 182520 : do mu=1,nd2gxdtfac
1294 [ - + ]: 84240 : if(cplex_d2gxdt(mu)==2)then
1295 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor)
1296 : : else
1297 [ + + ]: 252720 : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1298 : : end if
1299 : 84240 : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
1300 : 84240 : d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(1)
1301 : : #if !defined HAVE_OPENMP
1302 : 84240 : d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
1303 : 84240 : d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
1304 : : #endif
1305 [ + - ]: 168480 : if (cplex_==2) then
1306 : 84240 : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(2)
1307 : 84240 : d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
1308 : : #if !defined HAVE_OPENMP
1309 : 84240 : d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
1310 : 84240 : d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
1311 : : #endif
1312 : : end if
1313 : : end do
1314 : : end do
1315 : : #if defined HAVE_OPENMP
1316 : : if(jlmn<nlmn) then
1317 : : do ilmn=jlmn+1,nlmn
1318 : : i0lmn=ilmn*(ilmn-1)/2
1319 : : ijlmn=i0lmn+jlmn
1320 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
1321 : : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
1322 : : do mu=1,nd2gxdtfac
1323 : : if(cplex_d2gxdt(mu)==2)then
1324 : : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor)
1325 : : else
1326 : : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1327 : : end if
1328 : : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
1329 : : d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
1330 : : if (cplex_==2) then
1331 : : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
1332 : : d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
1333 : : end if
1334 : : end do
1335 : : end do
1336 : : end if
1337 : : #endif
1338 : : end do
1339 : : !$OMP END DO
1340 : : end do
1341 : : end do
1342 [ - + ]: 270 : ABI_FREE(gxfj)
1343 : : !$OMP END PARALLEL
1344 : : end if !nspinortot
1345 : : end if !complex
1346 : :
1347 : : ! === Off-diagonal term(s) (up-down, down-up)
1348 : :
1349 : : ! --- No parallelization over spinors ---
1350 [ + + + - ]: 4173 : if (nspinortot==2.and.nspinor==nspinortot) then
1351 [ - + ]: 270 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
1352 [ - + ]: 270 : ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!")
1353 : : !$OMP PARALLEL &
1354 : : !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), &
1355 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfi,gxfj,ilmn,i0lmn,ijlmn)
1356 [ - + + - : 1080 : ABI_MALLOC(gxfj,(cplex,nd2gxdtfac))
- + + - +
- + - - +
- + ]
1357 [ + + ]: 810 : do ispinor=1,nspinor
1358 : 540 : jspinor=3-ispinor
1359 [ + + ]: 1890 : do ia=1,nincat
1360 : 1080 : index_enl=atindx1(iatm+ia)
1361 : : !$OMP DO
1362 [ + + ]: 15660 : do jlmn=1,nlmn
1363 : 14040 : j0lmn=jlmn*(jlmn-1)/2
1364 : 14040 : jjlmn=j0lmn+jlmn
1365 [ + + ]: 42120 : enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor)
1366 [ + + ]: 28080 : do mu=1,nd2gxdtfac
1367 [ - + ]: 14040 : if(cplex_d2gxdt(mu)==2)then
1368 : 0 : cplex_ = 2 ;
1369 : 0 : gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,jlmn,ia,ispinor)
1370 : 0 : gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,jspinor)
1371 : : else
1372 : 14040 : cplex_ = cplex ;
1373 [ + + ]: 42120 : gxfi(1:cplex) =d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1374 [ + + ]: 42120 : gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,jspinor)
1375 : : end if
1376 : 14040 : d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
1377 : 14040 : d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
1378 [ + - ]: 28080 : if (cplex_==2) then
1379 : 14040 : d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
1380 : 14040 : d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
1381 : : end if
1382 : : end do
1383 [ + + ]: 99360 : do ilmn=1,jlmn-1
1384 : 84240 : ijlmn=j0lmn+ilmn
1385 [ + + ]: 252720 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
1386 [ + + ]: 182520 : do mu=1,nd2gxdtfac
1387 [ - + ]: 84240 : if(cplex_d2gxdt(mu)==2)then
1388 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor)
1389 : : else
1390 [ + + ]: 252720 : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1391 : : end if
1392 : 84240 : d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1)
1393 : 84240 : d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1)
1394 : : #if !defined HAVE_OPENMP
1395 : 84240 : d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu)
1396 : 84240 : d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu)
1397 : : #endif
1398 [ + - ]: 168480 : if (cplex_==2) then
1399 : 84240 : d2gxdtfac_(1,mu,jlmn,ia,jspinor)=d2gxdtfac_(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2)
1400 : 84240 : d2gxdtfac_(2,mu,jlmn,ia,jspinor)=d2gxdtfac_(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2)
1401 : : #if !defined HAVE_OPENMP
1402 : 84240 : d2gxdtfac_(1,mu,ilmn,ia,ispinor)=d2gxdtfac_(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu)
1403 : 84240 : d2gxdtfac_(2,mu,ilmn,ia,ispinor)=d2gxdtfac_(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu)
1404 : : #endif
1405 : : end if
1406 : : end do !mu
1407 : : end do !ilmn
1408 : : #if defined HAVE_OPENMP
1409 : : if(jlmn<nlmn) then
1410 : : do ilmn=jlmn+1,nlmn
1411 : : i0lmn=ilmn*(ilmn-1)/2
1412 : : ijlmn=i0lmn+jlmn
1413 : : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
1414 : : do mu=1,nd2gxdtfac
1415 : : if(cplex_d2gxdt(mu)==2)then
1416 : : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,jspinor)
1417 : : else
1418 : : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,jspinor)
1419 : : end if
1420 : : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1)
1421 : : d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1)
1422 : : if (cplex_==2) then
1423 : : d2gxdtfac_(1,mu,jlmn,ia,ispinor)=d2gxdtfac_(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2)
1424 : : d2gxdtfac_(2,mu,jlmn,ia,ispinor)=d2gxdtfac_(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2)
1425 : : end if
1426 : : end do !mu
1427 : : end do !ilmn
1428 : : end if
1429 : : #endif
1430 : : end do !jmln
1431 : : !$OMP END DO
1432 : : end do !ia
1433 : : end do !ispinor
1434 [ - + ]: 270 : ABI_FREE(gxfj)
1435 : : !$OMP END PARALLEL
1436 : :
1437 : : ! --- Parallelization over spinors ---
1438 [ - + - - ]: 3903 : else if (nspinortot==2.and.nspinor/=nspinortot) then
1439 [ # # ]: 0 : ABI_CHECK(cplex_enl==2,"BUG: opernlc_ylm: invalid cplex_enl/=2!")
1440 [ # # ]: 0 : ABI_CHECK(cplex_fac==2,"BUG: opernlc_ylm: invalid cplex_fac/=2!")
1441 [ # # # # : 0 : ABI_MALLOC(d2gxdtfac_offdiag,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinortot))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
1442 : : !$OMP PARALLEL &
1443 : : !$OMP PRIVATE(ia,index_enl), &
1444 : : !$OMP PRIVATE(jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,mu,gxfi)
1445 : : !$OMP WORKSHARE
1446 [ # # # # : 0 : d2gxdtfac_offdiag(:,:,:,:,:)=zero
# # # # #
# ]
1447 : : !$OMP END WORKSHARE
1448 : : !$OMP SINGLE
1449 : 0 : ispinor_index=mpi_enreg%me_spinor+1
1450 : 0 : jspinor_index=3-ispinor_index
1451 [ # # ]: 0 : if (ispinor_index==1) then
1452 : : ijspin=3;jispin=4
1453 : : else
1454 : 0 : ijspin=4;jispin=3
1455 : : end if
1456 : : !$OMP END SINGLE
1457 [ # # ]: 0 : do ia=1,nincat
1458 : 0 : index_enl=atindx1(iatm+ia)
1459 : : !$OMP DO
1460 [ # # ]: 0 : do jlmn=1,nlmn
1461 : 0 : j0lmn=jlmn*(jlmn-1)/2
1462 [ # # ]: 0 : do ilmn=1,nlmn
1463 : 0 : i0lmn=ilmn*(ilmn-1)/2
1464 [ # # ]: 0 : if (ilmn<=jlmn) then
1465 : 0 : ijlmn=j0lmn+ilmn
1466 : 0 : enl_(1)= enl_ptr(2*ijlmn-1,index_enl,ijspin)
1467 : 0 : enl_(2)=-enl_ptr(2*ijlmn ,index_enl,ijspin)
1468 : : else
1469 : 0 : jilmn=i0lmn+jlmn
1470 : 0 : enl_(1)= enl_ptr(2*jilmn-1,index_enl,jispin)
1471 : 0 : enl_(2)= enl_ptr(2*jilmn ,index_enl,jispin)
1472 : : end if
1473 [ # # ]: 0 : do mu=1,nd2gxdtfac
1474 [ # # ]: 0 : if(cplex_d2gxdt(mu)==2)then
1475 : 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1)
1476 : : else
1477 [ # # ]: 0 : cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1)
1478 : : end if
1479 : : d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
1480 : 0 : & d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(1)
1481 : : d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
1482 : 0 : & d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(2)*gxfi(1)
1483 [ # # ]: 0 : if (cplex_==2) then
1484 : : d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)= &
1485 : 0 : & d2gxdtfac_offdiag(1,mu,jlmn,ia,jspinor_index)-enl_(2)*gxfi(2)
1486 : : d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)= &
1487 : 0 : & d2gxdtfac_offdiag(2,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(2)
1488 : : end if
1489 : : end do
1490 : : end do !ilmn
1491 : : end do !jlmn
1492 : : !$OMP END DO
1493 : : end do !iat
1494 : : !$OMP SINGLE
1495 : 0 : call xmpi_sum(d2gxdtfac_offdiag,mpi_enreg%comm_spinor,ierr)
1496 : : !$OMP END SINGLE
1497 : : !$OMP WORKSHARE
1498 [ # # # # : 0 : d2gxdtfac_(:,:,:,:,1)=d2gxdtfac_(:,:,:,:,1)+d2gxdtfac_offdiag(:,:,:,:,ispinor_index)
# # # # ]
1499 : : !$OMP END WORKSHARE
1500 : : !$OMP END PARALLEL
1501 [ # # ]: 0 : ABI_FREE(d2gxdtfac_offdiag)
1502 : : end if !nspinortot
1503 : :
1504 : : end if ! pawopt & optder
1505 : :
1506 : : !End of loop when a exp(-iqR) phase is present
1507 : : !------------------------------------------- ------------------------
1508 : :
1509 : : !When iphase=1, gxfac and gxfac_ point to the same memory space
1510 : : !When iphase=2, we add i.gxfac_ to gxfac
1511 [ + + ]: 73991102 : if (iphase==2) then
1512 : : !$OMP PARALLEL PRIVATE(ispinor,ia,ilmn,mu)
1513 : : !$OMP DO COLLAPSE(3)
1514 [ + + ]: 10656280 : do ispinor=1,nspinor
1515 [ + + ]: 16057012 : do ia=1,nincat
1516 [ + + ]: 54328916 : do ilmn=1,nlmn
1517 : 43600044 : gxfac(1,ilmn,ia,ispinor)=gxfac(1,ilmn,ia,ispinor)-gxfac_(2,ilmn,ia,ispinor)
1518 : 49000776 : gxfac(2,ilmn,ia,ispinor)=gxfac(2,ilmn,ia,ispinor)+gxfac_(1,ilmn,ia,ispinor)
1519 : : end do
1520 : : end do
1521 : : end do
1522 : : !$OMP SINGLE
1523 [ - + ]: 5328140 : ABI_FREE(gxfac_)
1524 : : !$OMP END SINGLE
1525 [ + + ]: 5328140 : if (optder>=1) then
1526 : : !$OMP DO COLLAPSE(4)
1527 [ + + ]: 221688 : do ispinor=1,nspinor
1528 [ + + ]: 376740 : do ia=1,nincat
1529 [ + + ]: 1463940 : do ilmn=1,nlmn
1530 [ + + ]: 2551140 : do mu=1,ndgxdtfac
1531 : 1198044 : dgxdtfac(1,mu,ilmn,ia,ispinor)=dgxdtfac(1,mu,ilmn,ia,ispinor)-dgxdtfac_(2,mu,ilmn,ia,ispinor)
1532 : 2396088 : dgxdtfac(2,mu,ilmn,ia,ispinor)=dgxdtfac(2,mu,ilmn,ia,ispinor)+dgxdtfac_(1,mu,ilmn,ia,ispinor)
1533 : : end do
1534 : : end do
1535 : : end do
1536 : : end do
1537 : : !$OMP SINGLE
1538 [ - + ]: 110844 : ABI_FREE(dgxdtfac_)
1539 : : !$OMP END SINGLE
1540 : : end if
1541 [ - + ]: 5328140 : if (optder>=2) then
1542 : : !$OMP DO COLLAPSE(4)
1543 [ # # ]: 0 : do ispinor=1,nspinor
1544 [ # # ]: 0 : do ia=1,nincat
1545 [ # # ]: 0 : do ilmn=1,nlmn
1546 [ # # ]: 0 : do mu=1,nd2gxdtfac
1547 : 0 : d2gxdtfac(1,mu,ilmn,ia,ispinor)=d2gxdtfac(1,mu,ilmn,ia,ispinor)-d2gxdtfac_(2,mu,ilmn,ia,ispinor)
1548 : 0 : d2gxdtfac(2,mu,ilmn,ia,ispinor)=d2gxdtfac(2,mu,ilmn,ia,ispinor)+d2gxdtfac_(1,mu,ilmn,ia,ispinor)
1549 : : end do
1550 : : end do
1551 : : end do
1552 : : end do
1553 : : !$OMP SINGLE
1554 [ # # ]: 0 : ABI_FREE(d2gxdtfac_)
1555 : : !$OMP END SINGLE
1556 : : end if
1557 : : !$OMP END PARALLEL
1558 : : end if
1559 : :
1560 : : !End loop over real/imaginary part of the exp(-iqR) phase
1561 : : end do
1562 : :
1563 : :
1564 : : !Accumulate gxfac related to overlap (Sij) (PAW)
1565 : : !------------------------------------------- ------------------------
1566 [ + + ]: 35485077 : if (paw_opt==3.or.paw_opt==4) then ! Use Sij, overlap contribution
1567 : : !$OMP PARALLEL &
1568 : : !$OMP PRIVATE(ispinor,ia,jlmn,i0lmn,j0lmn,jjlmn,jlm,sijr,ilmn,ilm,ijlmn,gxi,gxj)
1569 : : !$OMP WORKSHARE
1570 [ + + + + : 946070713 : gxfac_sij(1:cplex,1:nlmn,1:nincat,1:nspinor)=zero
+ + + + ]
1571 : : !$OMP END WORKSHARE
1572 : : !$OMP DO COLLAPSE(3)
1573 [ + + ]: 37229480 : do ispinor=1,nspinor
1574 [ + + ]: 64962807 : do ia=1,nincat
1575 [ + + ]: 351831448 : do jlmn=1,nlmn
1576 : 304594418 : j0lmn=jlmn*(jlmn-1)/2
1577 : 304594418 : jjlmn=j0lmn+jlmn
1578 : 304594418 : jlm=indlmn(4,jlmn)
1579 : 304594418 : sijr=sij(jjlmn)
1580 [ + + ]: 881107906 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
1581 [ + + ]: 881107906 : gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxj(1:cplex)
1582 [ + + ]: 2198630912 : do ilmn=1,jlmn-1
1583 : 1866303167 : ilm=indlmn(4,ilmn)
1584 : : !if (ilm==jlm) then
1585 : 1866303167 : ijlmn=j0lmn+ilmn
1586 : 1866303167 : sijr=sij(ijlmn)
1587 [ + + ]: 5463783963 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
1588 [ + + ]: 5463783963 : gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex)
1589 : : #if !defined HAVE_OPENMP
1590 [ + + ]: 5768378381 : gxfac_sij(1:cplex,ilmn,ia,ispinor)=gxfac_sij(1:cplex,ilmn,ia,ispinor)+sijr*gxj(1:cplex)
1591 : : #endif
1592 : : !end if
1593 : : end do
1594 : : #if defined HAVE_OPENMP
1595 : : if(jlmn<nlmn) then
1596 : : do ilmn=jlmn+1,nlmn
1597 : : ilm=indlmn(4,ilmn)
1598 : : !if (ilm==jlm) then
1599 : : i0lmn=ilmn*(ilmn-1)/2
1600 : : ijlmn=i0lmn+jlmn
1601 : : sijr=sij(ijlmn)
1602 : : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
1603 : : gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex)
1604 : : !end if
1605 : : end do
1606 : : end if
1607 : : #endif
1608 : : end do
1609 : : end do
1610 : : end do
1611 : : !$OMP END DO
1612 : : !$OMP END PARALLEL
1613 : : end if
1614 : :
1615 : : !Accumulate dgxdtfac related to overlap (Sij) (PAW)
1616 : : !-------------------------------------------------------------------
1617 [ + + + + ]: 35485077 : if (optder>=1.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution
1618 : : !$OMP PARALLEL &
1619 : : !$OMP PRIVATE(ispinor,ia), &
1620 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,sijr,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
1621 [ - + + - : 7008692 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
- + + - +
- + - - +
- + ]
1622 : : !$OMP WORKSHARE
1623 [ + + + + : 80109388 : dgxdtfac_sij(1:cplex,1:ndgxdtfac,1:nlmn,1:nincat,1:nspinor)=zero
+ + + + +
+ ]
1624 : : !$OMP END WORKSHARE
1625 [ + + ]: 3593320 : do ispinor=1,nspinor
1626 [ + + ]: 5595978 : do ia=1,nincat
1627 : : !$OMP DO
1628 [ + + ]: 21380639 : do jlmn=1,nlmn
1629 : 17536834 : j0lmn=jlmn*(jlmn-1)/2
1630 : 17536834 : jjlmn=j0lmn+jlmn
1631 : 17536834 : sijr=sij(jjlmn)
1632 [ + + ]: 36529026 : do mu=1,ndgxdtfac
1633 [ + + ]: 56976576 : gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
1634 [ + + ]: 74513410 : dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1635 : : end do
1636 [ + + ]: 94742744 : do ilmn=1,jlmn-1
1637 : 75203252 : ijlmn=j0lmn+ilmn
1638 : 75203252 : sijr=sij(ijlmn)
1639 [ + + ]: 173936990 : do mu=1,ndgxdtfac
1640 [ + + ]: 243590712 : gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
1641 [ + + ]: 243590712 : dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1642 : : #if !defined HAVE_OPENMP
1643 [ + + ]: 318793964 : dgxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1644 : : #endif
1645 : : end do
1646 : : end do
1647 : : #if defined HAVE_OPENMP
1648 : : if(jlmn<nlmn) then
1649 : : do ilmn=jlmn+1,nlmn
1650 : : i0lmn=ilmn*(ilmn-1)/2
1651 : : ijlmn=i0lmn+jlmn
1652 : : sijr=sij(ijlmn)
1653 : : do mu=1,ndgxdtfac
1654 : : gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
1655 : : dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1656 : : end do
1657 : : end do
1658 : : end if
1659 : : #endif
1660 : : end do
1661 : : !$OMP END DO
1662 : : end do
1663 : : end do
1664 [ - + ]: 1752173 : ABI_FREE(gxfj)
1665 : : !$OMP END PARALLEL
1666 : : end if
1667 : :
1668 : : !Accumulate d2gxdtfac related to overlap (Sij) (PAW)
1669 : : !-------------------------------------------------------------------
1670 [ + + + + ]: 35485077 : if (optder==2.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution
1671 : : !$OMP PARALLEL &
1672 : : !$OMP PRIVATE(ispinor,ia), &
1673 : : !$OMP PRIVATE(jlmn,j0lmn,jjlmn,sijr,mu,gxfj,ilmn,i0lmn,ijlmn,gxfi)
1674 [ - + + - : 131892 : ABI_MALLOC(gxfj,(cplex,nd2gxdtfac))
- + + - +
- + - - +
- + ]
1675 : : !$OMP WORKSHARE
1676 [ + + + + : 1207254 : d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,1:nlmn,1:nincat,1:nspinor)=zero
+ + + + +
+ ]
1677 : : !$OMP END WORKSHARE
1678 [ + + ]: 66216 : do ispinor=1,nspinor
1679 [ + + ]: 100062 : do ia=1,nincat
1680 : : !$OMP DO
1681 [ + + ]: 343887 : do jlmn=1,nlmn
1682 : 276798 : j0lmn=jlmn*(jlmn-1)/2
1683 : 276798 : jjlmn=j0lmn+jlmn
1684 : 276798 : sijr=sij(jjlmn)
1685 [ + + ]: 553596 : do mu=1,nd2gxdtfac
1686 [ + + ]: 830394 : gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor)
1687 : : d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)= &
1688 [ + + ]: 1107192 : & d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1689 : : end do
1690 [ + + ]: 1318632 : do ilmn=1,jlmn-1
1691 : 1007988 : ijlmn=j0lmn+ilmn
1692 : 1007988 : sijr=sij(ijlmn)
1693 [ + + ]: 2292774 : do mu=1,nd2gxdtfac
1694 [ + + ]: 3023964 : gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1695 : : d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)= &
1696 [ + + ]: 3023964 : & d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1697 : : #if !defined HAVE_OPENMP
1698 : : d2gxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)= &
1699 [ + + ]: 4031952 : & d2gxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)+sijr*gxfj(1:cplex,mu)
1700 : : #endif
1701 : : end do
1702 : : end do
1703 : : #if defined HAVE_OPENMP
1704 : : if(jlmn<nlmn) then
1705 : : do ilmn=jlmn+1,nlmn
1706 : : i0lmn=ilmn*(ilmn-1)/2
1707 : : ijlmn=i0lmn+jlmn
1708 : : sijr=sij(ijlmn)
1709 : : do mu=1,nd2gxdtfac
1710 : : gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor)
1711 : : d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)= &
1712 : : & d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex)
1713 : : end do
1714 : : end do
1715 : : end if
1716 : : #endif
1717 : : end do
1718 : : !$OMP END DO
1719 : : end do
1720 : : end do
1721 [ - + ]: 32973 : ABI_FREE(gxfj)
1722 : : !$OMP END PARALLEL
1723 : : end if
1724 : :
1725 [ + + - + : 70970154 : end subroutine opernlc_ylm
- + - + -
+ ]
1726 : : !!***
1727 : :
1728 : : ! ---------------------------------------------------------------------------------------
1729 : :
1730 : : !!****f* m_opernlc_ylm/ls_ylm
1731 : : !! NAME
1732 : : !! ls_ylm
1733 : : !!
1734 : : !! FUNCTION
1735 : : !! Compute L.S operator matrix elements in real spherical harmonics basis.
1736 : : !! Upper triangle only (ilm<=jlm), packed as klm=jlm*(jlm-1)/2+ilm.
1737 : : !! ls_ylm(1,:,:)=Re, ls_ylm(2,:,:)=Im; ispin=1: up-up, ispin=2: up-dn.
1738 : : !! Adapted from m_paw_sphharm; tso debug blocks removed.
1739 : : !!
1740 : : !! SOURCE
1741 : :
1742 : 6570 : subroutine ls_ylm(ls_mat, lmax)
1743 : :
1744 : : !Arguments ---------------------------------------------
1745 : : integer, intent(in) :: lmax
1746 : : real(dp), allocatable, intent(inout) :: ls_mat(:,:,:)
1747 : :
1748 : : !Local variables ---------------------------------------
1749 : : integer :: im, jm, jlm, is, ll, lm0, mm
1750 : : real(dp), parameter :: isq2 = one/sqrt2
1751 : 6570 : complex(dp), allocatable :: U(:,:), LS(:,:,:), W(:,:)
1752 : : ! *************************************************************************
1753 : :
1754 [ + + + + : 1793610 : ls_mat = zero
+ + ]
1755 [ - + ]: 6570 : if (lmax <= 0) return
1756 : :
1757 [ + + ]: 19710 : do ll = 1, lmax
1758 : 13140 : lm0 = ll**2
1759 [ - + + - : 52560 : ABI_MALLOC(U, (2*ll+1, 2*ll+1))
- + + - +
- - + - +
- + ]
1760 [ - + + - : 65700 : ABI_MALLOC(LS, (2*ll+1, 2*ll+1, 2))
- + + - +
- + - - +
- + ]
1761 [ - + + - : 39420 : ABI_MALLOC(W, (2*ll+1, 2*ll+1))
- + + - -
+ - + ]
1762 [ + + + + : 867240 : U = czero; LS = czero
+ + + + +
+ ]
1763 : :
1764 : : ! Build U (real->complex Ylm transform) and LS (L.S in complex Ylm basis) in one pass
1765 [ + + ]: 65700 : do im = 1, 2*ll+1
1766 : 52560 : mm = im-ll-1
1767 [ + + ]: 52560 : if (mm > 0) then
1768 : 19710 : U(im,im) = (-1)**mm * isq2; U(-mm+ll+1,im) = isq2
1769 [ + + ]: 32850 : else if (mm == 0) then
1770 : 13140 : U(im,im) = cone
1771 : : else
1772 : 19710 : U(im,im) = cmplx(zero, isq2, dp); U(-mm+ll+1,im) = cmplx(zero, -(-1)**(-mm)*isq2, dp)
1773 : : end if
1774 : 52560 : LS(im,im,1) = half*mm
1775 [ + + ]: 52560 : if (mm+1 <= ll) LS(im,im+1,2) = half*sqrt(real((ll-mm)*(ll+mm+1), dp))
1776 [ + + ]: 65700 : if (mm-1 >= -ll) LS(im-1,im,2) = half*sqrt(real((ll+mm)*(ll-mm+1), dp))
1777 : : end do
1778 : :
1779 : : ! Transform to real Ylm basis via W = U^H * LS * U, store upper triangle
1780 [ + + ]: 39420 : do is = 1, 2
1781 [ - + + - : 8330760 : W = matmul(conjg(transpose(U)), matmul(LS(:,:,is), U))
- + + - +
- - + - +
+ + + + +
+ + + + +
+ - + - +
- - + - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - - - -
- - + + +
+ + + + +
+ + + - +
- - - - -
+ - ]
1782 [ + + ]: 144540 : do jm = 1, 2*ll+1
1783 : 105120 : jlm = lm0+jm
1784 [ + + ]: 407340 : do im = 1, jm
1785 : 275940 : ls_mat(1, jlm*(jlm-1)/2+lm0+im, is) = real(W(im,jm), dp)
1786 : 381060 : ls_mat(2, jlm*(jlm-1)/2+lm0+im, is) = aimag(W(im,jm))
1787 : : end do
1788 : : end do
1789 : : end do
1790 : :
1791 [ - + ]: 13140 : ABI_FREE(U)
1792 [ - + ]: 13140 : ABI_FREE(LS)
1793 [ - + ]: 19710 : ABI_FREE(W)
1794 : : end do
1795 : :
1796 [ - + ]: 6570 : end subroutine ls_ylm
1797 : : !!***
1798 : :
1799 [ + - ]: 26280 : end module m_opernlc_ylm
1800 : : !!***
|
