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