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 34745067 : subroutine opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,&
127 34745067 : & dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,dimekbq,enl,&
128 34745067 : & gx,gxfac,gxfac_sij,iatm,indlmn,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,&
129 34745067 : & 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 34745067 : real(dp), allocatable :: ls_ylm_so(:,:,:)
161 : !arrays
162 69490134 : real(dp) :: enl_(2),gxfi(2),gxi(cplex),gxj(cplex)
163 34745067 : real(dp),allocatable :: d2gxdtfac_offdiag(:,:,:,:,:),dgxdtfac_offdiag(:,:,:,:,:)
164 34745067 : real(dp),allocatable :: gxfac_offdiag(:,:,:,:),gxfj(:,:)
165 34745067 : real(dp),pointer :: d2gxdtfac_(:,:,:,:,:),dgxdtfac_(:,:,:,:,:),gxfac_(:,:,:,:)
166 34745067 : 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 74818274 : do iphase=1,dimekbq
179 40073207 : if (paw_opt==3) cycle
180 37794923 : if (iphase==1) then
181 32466783 : 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 1688282659 : gxfac_=zero
193 270832407 : if (optder>=1) dgxdtfac_=zero
194 108112184 : if (optder>=2) d2gxdtfac_=zero
195 37794923 : enl_ptr => enl(:,:,:,iphase)
196 :
197 : !NC+SO: precompute L.S matrix once (reused by gxfac, dgxdtfac, d2gxdtfac blocks below)
198 37794923 : lmax_so = 0
199 37794923 : 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 37794923 : if (paw_opt==0) then
215 : !Enl is E(Kleinman-Bylander)
216 9424629 : ABI_CHECK(cplex_enl/=2,"BUG: invalid cplex_enl=2!")
217 9424629 : 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 18862398 : do ispinor=1,nspinor
227 30685080 : do ia=1,nincat
228 143849433 : do ilmn=1,nlmn
229 122588982 : if (indlmn(6,ilmn)==2) cycle ! NC+SO: SO projectors handled separately below
230 122444442 : ispinor_index=ispinor+shift
231 122444442 : iln=indlmn(5,ilmn)
232 122444442 : enl_(1)=enl_ptr(iln,itypat,ispinor_index)
233 378778152 : 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 37794923 : 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 37794923 : 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 28370294 : 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 51339370 : do ispinor=1,nspinor
301 85788165 : do ia=1,nincat
302 399020977 : do jlmn=1,nlmn
303 338902497 : ispinor_index=ispinor+shift
304 338902497 : index_enl=atindx1(iatm+ia)
305 338902497 : j0lmn=jlmn*(jlmn-1)/2
306 338902497 : jjlmn=j0lmn+jlmn
307 338902497 : enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
308 338902497 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
309 981948121 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
310 981948121 : gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxj(1:cplex)
311 2201635845 : do ilmn=1,jlmn-1
312 1828284553 : ijlmn=j0lmn+ilmn
313 1828284553 : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
314 1828284553 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
315 5343255543 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
316 5343255543 : gxfac_(1:cplex,jlmn,ia,ispinor)=gxfac_(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex)
317 : #if !defined HAVE_OPENMP
318 5682158040 : 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 2700609 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
342 :
343 2700609 : 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 1704430 : do ia=1,nincat
348 852215 : index_enl=atindx1(iatm+ia)
349 : !$OMP DO
350 8522150 : do jlmn=1,nlmn
351 6817720 : j0lmn=jlmn*(jlmn-1)/2
352 6817720 : jjlmn=j0lmn+jlmn
353 6817720 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
354 6817720 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
355 20453160 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,1)
356 6817720 : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxj(1)
357 6817720 : if (cplex==2) gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxj(2)
358 31531955 : do ilmn=1,jlmn-1
359 23862020 : ijlmn=j0lmn+ilmn
360 71586060 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
361 23862020 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
362 71586060 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,1)
363 23862020 : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(1)*gxi(1)
364 23862020 : gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)-enl_(2)*gxi(1)
365 : #if !defined HAVE_OPENMP
366 23862020 : gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1)+enl_(1)*gxj(1)
367 23862020 : gxfac_(2,ilmn,ia,1)=gxfac_(2,ilmn,ia,1)+enl_(2)*gxj(1)
368 : #endif
369 30679740 : if (cplex==2) then
370 23862020 : gxfac_(1,jlmn,ia,1)=gxfac_(1,jlmn,ia,1)+enl_(2)*gxi(2)
371 23862020 : gxfac_(2,jlmn,ia,1)=gxfac_(2,jlmn,ia,1)+enl_(1)*gxi(2)
372 : #if !defined HAVE_OPENMP
373 23862020 : gxfac_(1,ilmn,ia,1)=gxfac_(1,ilmn,ia,1)-enl_(2)*gxj(2)
374 23862020 : 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 5484934 : do ispinor=1,nspinor
408 3636540 : ispinor_index=ispinor+shift
409 10059798 : do ia=1,nincat
410 4574864 : index_enl=atindx1(iatm+ia)
411 : !$OMP DO
412 71345604 : do jlmn=1,nlmn
413 63134200 : j0lmn=jlmn*(jlmn-1)/2
414 63134200 : jjlmn=j0lmn+jlmn
415 63134200 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,ispinor_index)
416 63134200 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
417 189402600 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
418 63134200 : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxj(1)
419 63134200 : if (cplex==2) gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxj(2)
420 522411720 : do ilmn=1,jlmn-1
421 454702656 : ijlmn=j0lmn+ilmn
422 1364107968 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index)
423 454702656 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
424 1364107968 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
425 454702656 : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(1)*gxi(1)
426 454702656 : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)-enl_(2)*gxi(1)
427 : #if !defined HAVE_OPENMP
428 454702656 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)+enl_(1)*gxj(1)
429 454702656 : gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(2)*gxj(1)
430 : #endif
431 517836856 : if (cplex==2) then
432 454702656 : gxfac_(1,jlmn,ia,ispinor)=gxfac_(1,jlmn,ia,ispinor)+enl_(2)*gxi(2)
433 454702656 : gxfac_(2,jlmn,ia,ispinor)=gxfac_(2,jlmn,ia,ispinor)+enl_(1)*gxi(2)
434 : #if !defined HAVE_OPENMP
435 454702656 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)-enl_(2)*gxj(2)
436 454702656 : 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 28370294 : if (nspinortot==2.and.nspinor==nspinortot) then
469 1788146 : ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!")
470 1788146 : 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 5364438 : do ispinor=1,nspinortot
475 3576292 : jspinor=3-ispinor
476 9836214 : do ia=1,nincat
477 4471776 : index_enl=atindx1(iatm+ia)
478 : !$OMP DO
479 69326684 : do jlmn=1,nlmn
480 61278616 : j0lmn=jlmn*(jlmn-1)/2
481 61278616 : jjlmn=j0lmn+jlmn
482 183835848 : enl_(1:2)=enl_ptr(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor )
483 183835848 : gxi(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
484 61278616 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(1)*gxi(1)
485 61278616 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)-enl_(2)*gxi(1)
486 61278616 : if (cplex==2) then
487 61278616 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(2)*gxi(2)
488 61278616 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)+enl_(1)*gxi(2)
489 : end if
490 : #if !defined HAVE_OPENMP
491 183835848 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,jspinor)
492 : #endif
493 504680584 : do ilmn=1,jlmn-1
494 438930192 : ijlmn=j0lmn+ilmn
495 1316790576 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor)
496 1316790576 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
497 438930192 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(1)*gxi(1)
498 438930192 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)-enl_(2)*gxi(1)
499 : #if !defined HAVE_OPENMP
500 438930192 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)+enl_(1)*gxj(1)
501 438930192 : gxfac_(2,ilmn,ia,ispinor)=gxfac_(2,ilmn,ia,ispinor)+enl_(2)*gxj(1)
502 : #endif
503 500208808 : if (cplex==2) then
504 438930192 : gxfac_(1,jlmn,ia,jspinor)=gxfac_(1,jlmn,ia,jspinor)+enl_(2)*gxi(2)
505 438930192 : gxfac_(2,jlmn,ia,jspinor)=gxfac_(2,jlmn,ia,jspinor)+enl_(1)*gxi(2)
506 : #if !defined HAVE_OPENMP
507 438930192 : gxfac_(1,ilmn,ia,ispinor)=gxfac_(1,ilmn,ia,ispinor)-enl_(2)*gxj(2)
508 438930192 : 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 26582148 : 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 37794923 : 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 37794923 : 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 37794923 : 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 1434789 : 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 4967196 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
672 2483598 : do ispinor=1,nspinor
673 1241799 : ispinor_index=ispinor+shift
674 3922388 : do ia=1,nincat
675 1438790 : index_enl=atindx1(iatm+ia)
676 : !$OMP DO
677 15282695 : do jlmn=1,nlmn
678 12602106 : j0lmn=jlmn*(jlmn-1)/2
679 12602106 : jjlmn=j0lmn+jlmn
680 12602106 : enl_(1)=enl_ptr(jjlmn,index_enl,ispinor_index)
681 12602106 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
682 25833434 : do mu=1,ndgxdtfac
683 39693984 : gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
684 52296090 : dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu)
685 : end do
686 68879912 : do ilmn=1,jlmn-1
687 54839016 : ijlmn=j0lmn+ilmn
688 54839016 : enl_(1)=enl_ptr(ijlmn,index_enl,ispinor_index)
689 54839016 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
690 125150126 : do mu=1,ndgxdtfac
691 173127012 : gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
692 173127012 : 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 227966028 : 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 1241799 : ABI_FREE(gxfj)
717 : !$OMP END PARALLEL
718 :
719 : ! 2-Enl is complex ===== D^ss'_ij=D^s's_ji^*
720 : else
721 192990 : ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!")
722 :
723 192990 : 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 489792 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
728 244896 : do ia=1,nincat
729 122448 : index_enl=atindx1(iatm+ia)
730 : !$OMP DO
731 1224480 : do jlmn=1,nlmn
732 979584 : j0lmn=jlmn*(jlmn-1)/2
733 979584 : jjlmn=j0lmn+jlmn
734 979584 : enl_(1)=enl_ptr(2*jjlmn-1,index_enl,1)
735 979584 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn)
736 1959168 : do mu=1,ndgxdtfac
737 979584 : 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 2938752 : cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,1)
741 : end if
742 979584 : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu)
743 1959168 : if (cplex_==2) dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu)
744 : end do
745 4530576 : do ilmn=1,jlmn-1
746 3428544 : ijlmn=j0lmn+ilmn
747 10285632 : enl_(1:2)=enl_ptr(2*ijlmn-1:2*ijlmn,index_enl,1)
748 3428544 : if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn)
749 7836672 : do mu=1,ndgxdtfac
750 3428544 : if(cplex_dgxdt(mu)==2)then
751 0 : cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1)
752 : else
753 10285632 : cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1)
754 : end if
755 3428544 : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1)
756 3428544 : dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)-enl_(2)*gxfi(1)
757 : #if !defined HAVE_OPENMP
758 3428544 : dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu)
759 3428544 : dgxdtfac_(2,mu,ilmn,ia,1)=dgxdtfac_(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu)
760 : #endif
761 6857088 : if (cplex_==2) then
762 3428544 : dgxdtfac_(1,mu,jlmn,ia,1)=dgxdtfac_(1,mu,jlmn,ia,1)+enl_(2)*gxfi(2)
763 3428544 : dgxdtfac_(2,mu,jlmn,ia,1)=dgxdtfac_(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2)
764 : #if !defined HAVE_OPENMP
765 3428544 : dgxdtfac_(1,mu,ilmn,ia,1)=dgxdtfac_(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu)
766 3428544 : 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 122448 : 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 1434789 : 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 1364247 : 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 37794923 : 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 37794923 : 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 37794923 : 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 37794923 : 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 72539990 : 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 34745067 : 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 921420637 : gxfac_sij(1:cplex,1:nlmn,1:nincat,1:nspinor)=zero
1571 : !$OMP END WORKSHARE
1572 : !$OMP DO COLLAPSE(3)
1573 35985600 : do ispinor=1,nspinor
1574 62848279 : do ia=1,nincat
1575 342113032 : do jlmn=1,nlmn
1576 296388174 : j0lmn=jlmn*(jlmn-1)/2
1577 296388174 : jjlmn=j0lmn+jlmn
1578 296388174 : jlm=indlmn(4,jlmn)
1579 296388174 : sijr=sij(jjlmn)
1580 858572358 : gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor)
1581 858572358 : gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxj(1:cplex)
1582 2152670676 : do ilmn=1,jlmn-1
1583 1829419823 : ilm=indlmn(4,ilmn)
1584 : !if (ilm==jlm) then
1585 1829419823 : ijlmn=j0lmn+ilmn
1586 1829419823 : sijr=sij(ijlmn)
1587 5363603795 : gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor)
1588 5363603795 : gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex)
1589 : #if !defined HAVE_OPENMP
1590 5659991969 : 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 34745067 : 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 6915956 : ABI_MALLOC(gxfj,(cplex,ndgxdtfac))
1622 : !$OMP WORKSHARE
1623 79067548 : dgxdtfac_sij(1:cplex,1:ndgxdtfac,1:nlmn,1:nincat,1:nspinor)=zero
1624 : !$OMP END WORKSHARE
1625 3546952 : do ispinor=1,nspinor
1626 5526426 : do ia=1,nincat
1627 : !$OMP DO
1628 21091199 : do jlmn=1,nlmn
1629 17293762 : j0lmn=jlmn*(jlmn-1)/2
1630 17293762 : jjlmn=j0lmn+jlmn
1631 17293762 : sijr=sij(jjlmn)
1632 36042882 : do mu=1,ndgxdtfac
1633 56247360 : gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor)
1634 73541122 : 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 93251336 : do ilmn=1,jlmn-1
1637 73978100 : ijlmn=j0lmn+ilmn
1638 73978100 : sijr=sij(ijlmn)
1639 171243614 : do mu=1,ndgxdtfac
1640 239915256 : gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor)
1641 239915256 : 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 313893356 : 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 1728989 : ABI_FREE(gxfj)
1665 : !$OMP END PARALLEL
1666 : end if
1667 :
1668 : !Accumulate d2gxdtfac related to overlap (Sij) (PAW)
1669 : !-------------------------------------------------------------------
1670 34745067 : 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 69490134 : 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 : !!***
|
