TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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)
enl(cplex_enl*dimenl1,dimenl2,nspinortot**2)=
->Norm conserving : ==== when paw_opt=0 ====
(Real) Kleinman-Bylander energies (hartree)
dimenl1=lmnmax - dimenl2=ntypat
->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)
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=informations 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)
hermdij=optional logical argument governing whether Dij are to be applied explicitly Hermitian (true)
or just symmetry Dij=Dji (false,default)
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.
PARENTS
m_gemm_nonlop,nonlop_ylm
CHILDREN
xmpi_sum
SOURCE
89 #if defined HAVE_CONFIG_H 90 #include "config.h" 91 #endif 92 93 #include "abi_common.h" 94 95 subroutine opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,dgxdt,dgxdtfac,dgxdtfac_sij,& 96 & d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,enl,gx,gxfac,gxfac_sij,iatm,indlmn,itypat,& 97 & lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,nincat,nlmn,& 98 & nspinor,nspinortot,optder,paw_opt,sij,hermdij) 99 100 use defs_basis 101 use defs_abitypes 102 use m_errors 103 use m_profiling_abi 104 use m_xmpi 105 106 !This section has been created automatically by the script Abilint (TD). 107 !Do not modify the following lines by hand. 108 #undef ABI_FUNC 109 #define ABI_FUNC 'opernlc_ylm' 110 !End of the abilint section 111 112 implicit none 113 114 !Arguments ------------------------------------ 115 !scalars 116 integer,intent(in) :: cplex,cplex_enl,cplex_fac,dimenl1,dimenl2,iatm,itypat 117 integer,intent(in) :: natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,nincat,nspinor,nspinortot,optder,paw_opt 118 integer,intent(inout) :: nlmn 119 real(dp) :: lambda 120 logical,optional,intent(in) :: hermdij 121 type(MPI_type) , intent(in) :: mpi_enreg 122 !arrays 123 integer,intent(in) :: atindx1(natom),indlmn(6,nlmn),cplex_dgxdt(ndgxdt),cplex_d2gxdt(nd2gxdt) 124 real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor) 125 real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor) 126 real(dp),intent(in) :: enl(dimenl1,dimenl2,nspinortot**2) 127 real(dp),intent(inout) :: gx(cplex,nlmn,nincat,nspinor) 128 real(dp),intent(in) :: sij(((paw_opt+1)/3)*nlmn*(nlmn+1)/2) 129 real(dp),intent(out) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor) 130 real(dp),intent(out) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3)) 131 real(dp),intent(out) :: d2gxdtfac(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor) 132 real(dp),intent(out) :: d2gxdtfac_sij(cplex,nd2gxdtfac,nlmn,nincat,nspinor*(paw_opt/3)) 133 real(dp),intent(out) :: gxfac(cplex_fac,nlmn,nincat,nspinor) 134 real(dp),intent(out) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3)) 135 136 !Local variables------------------------------- 137 !Arrays 138 !scalars 139 integer :: cplex_,ia,ierr,ijlmn,ijspin,ilm,ilmn,i0lmn,iln,index_enl,ispinor,ispinor_index 140 integer :: j0lmn,jilmn,jispin,jjlmn,jlm,jlmn,jspinor,jspinor_index,mu,shift 141 real(dp) :: sijr 142 logical :: hermdij_ 143 !arrays 144 real(dp) :: enl_(2),gxfi(2),gxi(cplex),gxj(cplex) 145 real(dp),allocatable :: d2gxdtfac_(:,:,:,:,:),dgxdtfac_(:,:,:,:,:),gxfac_(:,:,:,:),gxfj(:,:) 146 147 ! ************************************************************************* 148 149 DBG_ENTER("COLL") 150 151 !Parallelization over spinors treatment 152 shift=0;if (mpi_enreg%paral_spinor==1) shift=mpi_enreg%me_spinor 153 154 ! hermdij_ .false. invokes original coding, with Dij=Dji in complex case 155 ! hermdij_ .true. invokes Dij = Dji*, needed for magnetic field cases 156 hermdij_=.false.;if(present(hermdij)) hermdij_=hermdij 157 158 !Accumulate gxfac related to non-local operator (Norm-conserving) 159 !------------------------------------------------------------------- 160 if (paw_opt==0) then ! Enl is E(Kleinman-Bylander) 161 ABI_CHECK(cplex_enl/=2,"BUG: invalid cplex_enl=2!") 162 ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!") 163 !$OMP PARALLEL & 164 !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_) 165 !$OMP DO COLLAPSE(3) 166 do ispinor=1,nspinor 167 do ia=1,nincat 168 do ilmn=1,nlmn 169 ispinor_index=ispinor+shift 170 iln=indlmn(5,ilmn) 171 enl_(1)=enl(iln,itypat,ispinor_index) 172 gxfac(1:cplex,ilmn,ia,ispinor)=enl_(1)*gx(1:cplex,ilmn,ia,ispinor) 173 end do 174 end do 175 end do 176 !$OMP END DO 177 !$OMP END PARALLEL 178 end if 179 180 !Accumulate gxfac related to nonlocal operator (PAW) 181 !------------------------------------------------------------------- 182 if (paw_opt==1.or.paw_opt==2.or.paw_opt==4) then ! Enl is psp strength Dij 183 gxfac(1:cplex_fac,1:nlmn,1:nincat,1:nspinor)=zero ! or (Dij-lambda.Sij) 184 ! === Diagonal term(s) (up-up, down-down) 185 ! 1-Enl is real 186 if (cplex_enl==1) then 187 188 !$OMP PARALLEL & 189 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), & 190 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,ijlmn,gxi) 191 !$OMP DO COLLAPSE(3) 192 do ispinor=1,nspinor 193 do ia=1,nincat 194 do jlmn=1,nlmn 195 ispinor_index=ispinor+shift 196 index_enl=atindx1(iatm+ia) 197 j0lmn=jlmn*(jlmn-1)/2 198 jjlmn=j0lmn+jlmn 199 enl_(1)=enl(jjlmn,index_enl,ispinor_index) 200 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 201 gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor) 202 gxfac(1:cplex,jlmn,ia,ispinor)=gxfac(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxj(1:cplex) 203 do ilmn=1,jlmn-1 204 ijlmn=j0lmn+ilmn 205 enl_(1)=enl(ijlmn,index_enl,ispinor_index) 206 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 207 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 208 gxfac(1:cplex,jlmn,ia,ispinor)=gxfac(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex) 209 #if !defined HAVE_OPENMP 210 gxfac(1:cplex,ilmn,ia,ispinor)=gxfac(1:cplex,ilmn,ia,ispinor)+enl_(1)*gxj(1:cplex) 211 #endif 212 end do 213 #if defined HAVE_OPENMP 214 if(jlmn<nlmn) then 215 do ilmn=jlmn+1,nlmn 216 i0lmn=(ilmn*(ilmn-1)/2) 217 ijlmn=i0lmn+jlmn 218 enl_(1)=enl(ijlmn,index_enl,ispinor_index) 219 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 220 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 221 gxfac(1:cplex,jlmn,ia,ispinor)=gxfac(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxi(1:cplex) 222 end do 223 end if 224 #endif 225 end do 226 end do 227 end do 228 !$OMP END DO 229 !$OMP END PARALLEL 230 231 ! 2-Enl is complex 232 else 233 ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!") 234 if (nspinortot==1) then !===== when nspinor=1, D_ij=D_ji except in hermdij case 235 236 !$OMP PARALLEL & 237 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,ijlmn,gxi) 238 do ia=1,nincat 239 index_enl=atindx1(iatm+ia) 240 !$OMP DO 241 do jlmn=1,nlmn 242 j0lmn=jlmn*(jlmn-1)/2 243 jjlmn=j0lmn+jlmn 244 enl_(1:2)=enl(2*jjlmn-1:2*jjlmn,index_enl,1) 245 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 246 gxj(1:cplex)=gx(1:cplex,jlmn,ia,1) 247 if(hermdij_) then 248 gxfac(1:cplex,jlmn,ia,1)=gxfac(1:cplex,jlmn,ia,1)+enl_(1)*gxj(1:cplex) 249 else 250 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)+enl_(1)*gxj(1) 251 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)+enl_(2)*gxj(1) 252 if (cplex==2) then 253 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)-enl_(2)*gxj(2) 254 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)+enl_(1)*gxj(2) 255 end if 256 end if 257 do ilmn=1,jlmn-1 258 ijlmn=j0lmn+ilmn 259 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,1) 260 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 261 gxi(1:cplex)=gx(1:cplex,ilmn,ia,1) 262 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)+enl_(1)*gxi(1) 263 if(hermdij_) then 264 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)-enl_(2)*gxi(1) 265 else 266 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)+enl_(2)*gxi(1) 267 end if 268 #if !defined HAVE_OPENMP 269 gxfac(1,ilmn,ia,1)=gxfac(1,ilmn,ia,1)+enl_(1)*gxj(1) 270 gxfac(2,ilmn,ia,1)=gxfac(2,ilmn,ia,1)+enl_(2)*gxj(1) 271 #endif 272 if (cplex==2) then 273 if(hermdij_) then 274 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)+enl_(2)*gxi(2) 275 else 276 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)-enl_(2)*gxi(2) 277 end if 278 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)+enl_(1)*gxi(2) 279 #if !defined HAVE_OPENMP 280 gxfac(1,ilmn,ia,1)=gxfac(1,ilmn,ia,1)-enl_(2)*gxj(2) 281 gxfac(2,ilmn,ia,1)=gxfac(2,ilmn,ia,1)+enl_(1)*gxj(2) 282 #endif 283 end if 284 end do 285 #if defined HAVE_OPENMP 286 if(jlmn<nlmn) then 287 do ilmn=jlmn+1,nlmn 288 i0lmn=ilmn*(ilmn-1)/2 289 ijlmn=i0lmn+jlmn 290 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,1) 291 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 292 gxi(1:cplex)=gx(1:cplex,ilmn,ia,1) 293 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)+enl_(1)*gxi(1) 294 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)+enl_(2)*gxi(1) 295 if (cplex==2) then 296 gxfac(1,jlmn,ia,1)=gxfac(1,jlmn,ia,1)-enl_(2)*gxi(2) 297 gxfac(2,jlmn,ia,1)=gxfac(2,jlmn,ia,1)+enl_(1)*gxi(2) 298 end if 299 end do 300 end if 301 #endif 302 end do 303 !$OMP END DO 304 end do 305 !$OMP END PARALLEL 306 else !===== when nspinor=2, D_ij=D_ji^* 307 !$OMP PARALLEL & 308 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), & 309 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxj,ilmn,ijlmn,gxi) 310 do ispinor=1,nspinor 311 ispinor_index=ispinor+shift 312 do ia=1,nincat 313 index_enl=atindx1(iatm+ia) 314 !$OMP DO 315 do jlmn=1,nlmn 316 j0lmn=jlmn*(jlmn-1)/2 317 jjlmn=j0lmn+jlmn 318 enl_(1)=enl(2*jjlmn-1,index_enl,ispinor_index) ! enl_ii is real for up-up and dn-dn 319 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 320 gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor) 321 gxfac(1:cplex,jlmn,ia,ispinor)=gxfac(1:cplex,jlmn,ia,ispinor)+enl_(1)*gxj(1:cplex) 322 do ilmn=1,jlmn-1 323 ijlmn=j0lmn+ilmn 324 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index) 325 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 326 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 327 gxfac(1,jlmn,ia,ispinor)=gxfac(1,jlmn,ia,ispinor)+enl_(1)*gxi(1) 328 gxfac(2,jlmn,ia,ispinor)=gxfac(2,jlmn,ia,ispinor)-enl_(2)*gxi(1) 329 #if !defined HAVE_OPENMP 330 gxfac(1,ilmn,ia,ispinor)=gxfac(1,ilmn,ia,ispinor)+enl_(1)*gxj(1) 331 gxfac(2,ilmn,ia,ispinor)=gxfac(2,ilmn,ia,ispinor)+enl_(2)*gxj(1) 332 #endif 333 if (cplex==2) then 334 gxfac(1,jlmn,ia,ispinor)=gxfac(1,jlmn,ia,ispinor)+enl_(2)*gxi(2) 335 gxfac(2,jlmn,ia,ispinor)=gxfac(2,jlmn,ia,ispinor)+enl_(1)*gxi(2) 336 #if !defined HAVE_OPENMP 337 gxfac(1,ilmn,ia,ispinor)=gxfac(1,ilmn,ia,ispinor)-enl_(2)*gxj(2) 338 gxfac(2,ilmn,ia,ispinor)=gxfac(2,ilmn,ia,ispinor)+enl_(1)*gxj(2) 339 #endif 340 end if 341 end do 342 #if defined HAVE_OPENMP 343 if(jlmn<nlmn) then 344 do ilmn=jlmn+1,nlmn 345 i0lmn=ilmn*(ilmn-1)/2 346 ijlmn=i0lmn+jlmn 347 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index) 348 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 349 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 350 gxfac(1,jlmn,ia,ispinor)=gxfac(1,jlmn,ia,ispinor)+enl_(1)*gxi(1) 351 gxfac(2,jlmn,ia,ispinor)=gxfac(2,jlmn,ia,ispinor)+enl_(2)*gxi(1) 352 if (cplex==2) then 353 gxfac(1,jlmn,ia,ispinor)=gxfac(1,jlmn,ia,ispinor)-enl_(2)*gxi(2) 354 gxfac(2,jlmn,ia,ispinor)=gxfac(2,jlmn,ia,ispinor)+enl_(1)*gxi(2) 355 end if 356 end do 357 end if 358 #endif 359 end do 360 !$OMP END DO 361 end do 362 end do 363 !$OMP END PARALLEL 364 end if !nspinortot 365 end if !complex_enl 366 367 ! === Off-diagonal term(s) (up-down, down-up) 368 ! --- No parallelization over spinors --- 369 if (nspinortot==2.and.nspinor==nspinortot) then 370 ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!") 371 ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex)!") 372 !$OMP PARALLEL & 373 !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), & 374 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,gxi,gxj,ilmn,ijlmn) 375 do ispinor=1,nspinortot 376 jspinor=3-ispinor 377 do ia=1,nincat 378 index_enl=atindx1(iatm+ia) 379 !$OMP DO 380 do jlmn=1,nlmn 381 j0lmn=jlmn*(jlmn-1)/2 382 jjlmn=j0lmn+jlmn 383 enl_(1:2)=enl(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor ) 384 gxi(1:cplex)=gx(1:cplex,jlmn,ia,ispinor) 385 gxfac(1,jlmn,ia,jspinor)=gxfac(1,jlmn,ia,jspinor)+enl_(1)*gxi(1) 386 gxfac(2,jlmn,ia,jspinor)=gxfac(2,jlmn,ia,jspinor)-enl_(2)*gxi(1) 387 if (cplex==2) then 388 gxfac(1,jlmn,ia,jspinor)=gxfac(1,jlmn,ia,jspinor)+enl_(2)*gxi(2) 389 gxfac(2,jlmn,ia,jspinor)=gxfac(2,jlmn,ia,jspinor)+enl_(1)*gxi(2) 390 end if 391 #if !defined HAVE_OPENMP 392 gxj(1:cplex)=gx(1:cplex,jlmn,ia,jspinor) 393 #endif 394 do ilmn=1,jlmn-1 395 ijlmn=j0lmn+ilmn 396 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor) 397 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 398 gxfac(1,jlmn,ia,jspinor)=gxfac(1,jlmn,ia,jspinor)+enl_(1)*gxi(1) 399 gxfac(2,jlmn,ia,jspinor)=gxfac(2,jlmn,ia,jspinor)-enl_(2)*gxi(1) 400 #if !defined HAVE_OPENMP 401 gxfac(1,ilmn,ia,ispinor)=gxfac(1,ilmn,ia,ispinor)+enl_(1)*gxj(1) 402 gxfac(2,ilmn,ia,ispinor)=gxfac(2,ilmn,ia,ispinor)+enl_(2)*gxj(1) 403 #endif 404 if (cplex==2) then 405 gxfac(1,jlmn,ia,jspinor)=gxfac(1,jlmn,ia,jspinor)+enl_(2)*gxi(2) 406 gxfac(2,jlmn,ia,jspinor)=gxfac(2,jlmn,ia,jspinor)+enl_(1)*gxi(2) 407 #if !defined HAVE_OPENMP 408 gxfac(1,ilmn,ia,ispinor)=gxfac(1,ilmn,ia,ispinor)-enl_(2)*gxj(2) 409 gxfac(2,ilmn,ia,ispinor)=gxfac(2,ilmn,ia,ispinor)+enl_(1)*gxj(2) 410 #endif 411 end if 412 end do 413 #if defined HAVE_OPENMP 414 if(jlmn<nlmn) then 415 do ilmn=jlmn+1,nlmn 416 i0lmn=ilmn*(ilmn-1)/2 417 ijlmn=i0lmn+jlmn 418 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor) 419 gxi(1:cplex)=gx(1:cplex,ilmn,ia,jspinor) 420 gxfac(1,jlmn,ia,ispinor)=gxfac(1,jlmn,ia,ispinor)+enl_(1)*gxi(1) 421 gxfac(2,jlmn,ia,ispinor)=gxfac(2,jlmn,ia,ispinor)+enl_(2)*gxi(1) 422 if (cplex==2) then 423 gxfac(1,jlmn,ia,ispinor)=gxfac(1,jlmn,ia,ispinor)-enl_(2)*gxi(2) 424 gxfac(2,jlmn,ia,ispinor)=gxfac(2,jlmn,ia,ispinor)+enl_(1)*gxi(2) 425 end if 426 end do 427 end if 428 #endif 429 end do 430 !$OMP END DO 431 end do 432 end do 433 !$OMP END PARALLEL 434 ! --- Parallelization over spinors --- 435 else if (nspinortot==2.and.nspinor/=nspinortot) then 436 ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!") 437 ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!") 438 ABI_ALLOCATE(gxfac_,(cplex_fac,nlmn,nincat,nspinortot)) 439 gxfac_(:,:,:,:)=zero 440 ispinor_index=mpi_enreg%me_spinor+1 441 jspinor_index=3-ispinor_index 442 if (ispinor_index==1) then 443 ijspin=3;jispin=4 444 else 445 ijspin=4;jispin=3 446 end if 447 !$OMP PARALLEL & 448 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,gxi) 449 do ia=1,nincat 450 index_enl=atindx1(iatm+ia) 451 !$OMP DO 452 do jlmn=1,nlmn 453 j0lmn=jlmn*(jlmn-1)/2 454 do ilmn=1,nlmn 455 i0lmn=ilmn*(ilmn-1)/2 456 if (ilmn<=jlmn) then 457 ijlmn=j0lmn+ilmn 458 enl_(1)= enl(2*ijlmn-1,index_enl,ijspin) 459 enl_(2)=-enl(2*ijlmn ,index_enl,ijspin) 460 else 461 jilmn=i0lmn+jlmn 462 enl_(1:2)=enl(2*jilmn-1:2*jilmn,index_enl,jispin) 463 end if 464 gxi(1:cplex)=gx(1:cplex,ilmn,ia,1) 465 gxfac_(1,jlmn,ia,jspinor_index)=gxfac_(1,jlmn,ia,jspinor_index)+enl_(1)*gxi(1) 466 gxfac_(2,jlmn,ia,jspinor_index)=gxfac_(2,jlmn,ia,jspinor_index)+enl_(2)*gxi(1) 467 if (cplex==2) then 468 gxfac_(1,jlmn,ia,jspinor_index)=gxfac_(1,jlmn,ia,jspinor_index)-enl_(2)*gxi(2) 469 gxfac_(2,jlmn,ia,jspinor_index)=gxfac_(2,jlmn,ia,jspinor_index)+enl_(1)*gxi(2) 470 end if 471 end do !ilmn 472 end do !jlmn 473 !$OMP END DO 474 end do !iat 475 !$OMP END PARALLEL 476 call xmpi_sum(gxfac_,mpi_enreg%comm_spinor,ierr) 477 gxfac(:,:,:,1)=gxfac(:,:,:,1)+gxfac_(:,:,:,ispinor_index) 478 ABI_DEALLOCATE(gxfac_) 479 end if 480 481 end if !paw_opt 482 483 !Accumulate gxfac related to overlap (Sij) (PAW) 484 !------------------------------------------- ------------------------ 485 if (paw_opt==3.or.paw_opt==4) then ! Use Sij, overlap contribution 486 !$OMP PARALLEL & 487 !$OMP PRIVATE(ispinor,ia,jlmn,j0lmn,jjlmn,jlm,sijr,ilmn,ilm,ijlmn,gxi,gxj) 488 !$OMP WORKSHARE 489 gxfac_sij(1:cplex,1:nlmn,1:nincat,1:nspinor)=zero 490 !$OMP END WORKSHARE 491 !$OMP DO COLLAPSE(3) 492 do ispinor=1,nspinor 493 do ia=1,nincat 494 do jlmn=1,nlmn 495 j0lmn=jlmn*(jlmn-1)/2 496 jjlmn=j0lmn+jlmn 497 jlm=indlmn(4,jlmn) 498 sijr=sij(jjlmn) 499 gxj(1:cplex)=gx(1:cplex,jlmn,ia,ispinor) 500 gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxj(1:cplex) 501 do ilmn=1,jlmn-1 502 ilm=indlmn(4,ilmn) 503 !if (ilm==jlm) then 504 ijlmn=j0lmn+ilmn 505 sijr=sij(ijlmn) 506 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 507 gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex) 508 #if !defined HAVE_OPENMP 509 gxfac_sij(1:cplex,ilmn,ia,ispinor)=gxfac_sij(1:cplex,ilmn,ia,ispinor)+sijr*gxj(1:cplex) 510 #endif 511 !end if 512 end do 513 #if defined HAVE_OPENMP 514 if(jlmn<nlmn) then 515 do ilmn=jlmn+1,nlmn 516 ilm=indlmn(4,ilmn) 517 !if (ilm==jlm) then 518 i0lmn=ilmn*(ilmn-1)/2 519 ijlmn=i0lmn+jlmn 520 sijr=sij(ijlmn) 521 gxi(1:cplex)=gx(1:cplex,ilmn,ia,ispinor) 522 gxfac_sij(1:cplex,jlmn,ia,ispinor)=gxfac_sij(1:cplex,jlmn,ia,ispinor)+sijr*gxi(1:cplex) 523 !end if 524 end do 525 end if 526 #endif 527 end do 528 end do 529 end do 530 !$OMP END DO 531 !$OMP END PARALLEL 532 end if 533 534 !Accumulate dgxdtfac related to nonlocal operator (Norm-conserving) 535 !------------------------------------------------------------------- 536 if (optder>=1.and.paw_opt==0) then ! Enl is E(Kleinman-Bylander) 537 ABI_CHECK(cplex_enl==1,"BUG: invalid cplex_enl/=1!") 538 ABI_CHECK(cplex_fac==cplex,"BUG: invalid cplex_fac/=cplex!") 539 !$OMP PARALLEL & 540 !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_,mu) 541 do ispinor=1,nspinor 542 ispinor_index = ispinor + shift 543 do ia=1,nincat 544 !$OMP DO 545 do ilmn=1,nlmn 546 iln=indlmn(5,ilmn) 547 enl_(1)=enl(iln,itypat,ispinor_index) 548 do mu=1,ndgxdtfac 549 dgxdtfac(1:cplex,mu,ilmn,ia,ispinor)=enl_(1)*dgxdt(1:cplex,mu,ilmn,ia,ispinor) 550 end do 551 end do 552 !$OMP END DO 553 end do 554 end do 555 !$OMP END PARALLEL 556 end if 557 558 !Accumulate dgxdtfac related to nonlocal operator (PAW) 559 !------------------------------------------------------------------- 560 if (optder>=1.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then ! Enl is psp strength Dij 561 dgxdtfac(1:cplex,1:ndgxdtfac,1:nlmn,1:nincat,1:nspinor)=zero 562 ! === Diagonal term(s) (up-up, down-down) 563 ! 1-Enl is real 564 if (cplex_enl==1) then 565 !$OMP PARALLEL & 566 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi) 567 ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac)) 568 do ispinor=1,nspinor 569 ispinor_index=ispinor+shift 570 do ia=1,nincat 571 index_enl=atindx1(iatm+ia) 572 !$OMP DO 573 do jlmn=1,nlmn 574 j0lmn=jlmn*(jlmn-1)/2 575 jjlmn=j0lmn+jlmn 576 enl_(1)=enl(jjlmn,index_enl,ispinor_index) 577 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 578 do mu=1,ndgxdtfac 579 gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor) 580 dgxdtfac(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu) 581 end do 582 do ilmn=1,jlmn-1 583 ijlmn=j0lmn+ilmn 584 enl_(1)=enl(ijlmn,index_enl,ispinor_index) 585 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 586 do mu=1,ndgxdtfac 587 gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 588 dgxdtfac(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex) 589 #if !defined HAVE_OPENMP 590 dgxdtfac(1:cplex,mu,ilmn,ia,ispinor)=dgxdtfac(1:cplex,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu) 591 #endif 592 end do 593 end do 594 #if defined HAVE_OPENMP 595 if(jlmn<nlmn) then 596 do ilmn=jlmn+1,nlmn 597 i0lmn=ilmn*(ilmn-1)/2 598 ijlmn=i0lmn+jlmn 599 enl_(1)=enl(ijlmn,index_enl,ispinor_index) 600 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 601 do mu=1,ndgxdtfac 602 gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 603 dgxdtfac(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex) 604 end do 605 end do 606 end if 607 #endif 608 end do 609 !$OMP END DO 610 end do 611 end do 612 ABI_DEALLOCATE(gxfj) 613 !$OMP END PARALLEL 614 ! 2-Enl is complex 615 else 616 ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!") 617 if (nspinortot==1) then !===== when nspinor=1, D_ij=D_ji 618 !$OMP PARALLEL & 619 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi) 620 ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac)) 621 do ia=1,nincat 622 index_enl=atindx1(iatm+ia) 623 !$OMP DO 624 do jlmn=1,nlmn 625 j0lmn=jlmn*(jlmn-1)/2 626 jjlmn=j0lmn+jlmn 627 enl_(1:2)=enl(2*jjlmn-1:2*jjlmn,index_enl,1) 628 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 629 do mu=1,ndgxdtfac 630 if(cplex_dgxdt(mu)==2)then 631 cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,1) 632 else 633 cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,1) 634 end if 635 dgxdtfac(1,mu,jlmn,ia,1)=dgxdtfac(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu) 636 dgxdtfac(2,mu,jlmn,ia,1)=dgxdtfac(2,mu,jlmn,ia,1)+enl_(2)*gxfj(1,mu) 637 if (cplex_==2) then 638 dgxdtfac(1,mu,jlmn,ia,1)=dgxdtfac(1,mu,jlmn,ia,1)-enl_(2)*gxfj(2,mu) 639 dgxdtfac(2,mu,jlmn,ia,1)=dgxdtfac(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu) 640 end if 641 end do 642 do ilmn=1,jlmn-1 643 ijlmn=j0lmn+ilmn 644 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,1) 645 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 646 do mu=1,ndgxdtfac 647 if(cplex_dgxdt(mu)==2)then 648 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1) 649 else 650 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1) 651 end if 652 dgxdtfac(1,mu,jlmn,ia,1)=dgxdtfac(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1) 653 dgxdtfac(2,mu,jlmn,ia,1)=dgxdtfac(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1) 654 #if !defined HAVE_OPENMP 655 dgxdtfac(1,mu,ilmn,ia,1)=dgxdtfac(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu) 656 dgxdtfac(2,mu,ilmn,ia,1)=dgxdtfac(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu) 657 #endif 658 if (cplex_==2) then 659 dgxdtfac(1,mu,jlmn,ia,1)=dgxdtfac(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2) 660 dgxdtfac(2,mu,jlmn,ia,1)=dgxdtfac(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2) 661 #if !defined HAVE_OPENMP 662 dgxdtfac(1,mu,ilmn,ia,1)=dgxdtfac(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu) 663 dgxdtfac(2,mu,ilmn,ia,1)=dgxdtfac(2,mu,ilmn,ia,1)+enl_(1)*gxfj(2,mu) 664 #endif 665 end if 666 end do 667 end do 668 #if defined HAVE_OPENMP 669 if(jlmn<nlmn) then 670 do ilmn=jlmn+1,nlmn 671 i0lmn=ilmn*(ilmn-1)/2 672 ijlmn=i0lmn+jlmn 673 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,1) 674 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 675 do mu=1,ndgxdtfac 676 if(cplex_dgxdt(mu)==2)then 677 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1) 678 else 679 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1) 680 end if 681 dgxdtfac(1,mu,jlmn,ia,1)=dgxdtfac(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1) 682 dgxdtfac(2,mu,jlmn,ia,1)=dgxdtfac(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1) 683 if (cplex_==2) then 684 dgxdtfac(1,mu,jlmn,ia,1)=dgxdtfac(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2) 685 dgxdtfac(2,mu,jlmn,ia,1)=dgxdtfac(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2) 686 end if 687 end do 688 end do 689 end if 690 #endif 691 end do 692 !$OMP END DO 693 end do 694 ABI_DEALLOCATE(gxfj) 695 !$OMP END PARALLEL 696 else !===== when nspinor=2, D_ij=D_ji^* 697 !$OMP PARALLEL & 698 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), & 699 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi) 700 ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac)) 701 do ispinor=1,nspinor 702 ispinor_index = ispinor + shift 703 do ia=1,nincat 704 index_enl=atindx1(iatm+ia) 705 !$OMP DO 706 do jlmn=1,nlmn 707 j0lmn=jlmn*(jlmn-1)/2 708 jjlmn=j0lmn+jlmn 709 enl_(1)=enl(2*jjlmn-1,index_enl,ispinor_index) ! enl_ii is real for up-up and dn-dn 710 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 711 do mu=1,ndgxdtfac 712 if(cplex_dgxdt(mu)==2)then 713 cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,ispinor) 714 else 715 cplex_ = cplex ; gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor) 716 end if 717 dgxdtfac(1:cplex_,mu,jlmn,ia,ispinor)=dgxdtfac(1:cplex_,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex_,mu) 718 end do 719 do ilmn=1,jlmn-1 720 ijlmn=j0lmn+ilmn 721 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index) 722 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 723 do mu=1,ndgxdtfac 724 if(cplex_dgxdt(mu)==2)then 725 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor) 726 else 727 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 728 end if 729 dgxdtfac(1,mu,jlmn,ia,ispinor)=dgxdtfac(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1) 730 dgxdtfac(2,mu,jlmn,ia,ispinor)=dgxdtfac(2,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(1) 731 #if !defined HAVE_OPENMP 732 dgxdtfac(1,mu,ilmn,ia,ispinor)=dgxdtfac(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu) 733 dgxdtfac(2,mu,ilmn,ia,ispinor)=dgxdtfac(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu) 734 #endif 735 if (cplex_==2) then 736 dgxdtfac(1,mu,jlmn,ia,ispinor)=dgxdtfac(1,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(2) 737 dgxdtfac(2,mu,jlmn,ia,ispinor)=dgxdtfac(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2) 738 #if !defined HAVE_OPENMP 739 dgxdtfac(1,mu,ilmn,ia,ispinor)=dgxdtfac(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu) 740 dgxdtfac(2,mu,ilmn,ia,ispinor)=dgxdtfac(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu) 741 #endif 742 end if 743 end do 744 end do 745 #if defined HAVE_OPENMP 746 if(jlmn<nlmn) then 747 do ilmn=jlmn+1,nlmn 748 i0lmn=ilmn*(ilmn-1)/2 749 ijlmn=i0lmn+jlmn 750 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index) 751 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 752 do mu=1,ndgxdtfac 753 if(cplex_dgxdt(mu)==2)then 754 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor) 755 else 756 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 757 end if 758 dgxdtfac(1,mu,jlmn,ia,ispinor)=dgxdtfac(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1) 759 dgxdtfac(2,mu,jlmn,ia,ispinor)=dgxdtfac(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1) 760 if (cplex_==2) then 761 dgxdtfac(1,mu,jlmn,ia,ispinor)=dgxdtfac(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2) 762 dgxdtfac(2,mu,jlmn,ia,ispinor)=dgxdtfac(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2) 763 end if 764 end do 765 end do 766 end if 767 #endif 768 end do 769 !$OMP END DO 770 end do 771 end do 772 ABI_DEALLOCATE(gxfj) 773 !$OMP END PARALLEL 774 end if !nspinortot 775 end if !complex 776 ! === Off-diagonal term(s) (up-down, down-up) 777 ! --- No parallelization over spinors --- 778 if (nspinortot==2.and.nspinor==nspinortot) then 779 ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!") 780 ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!") 781 !$OMP PARALLEL & 782 !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), & 783 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfi,gxfj,ilmn,ijlmn) 784 ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac)) 785 do ispinor=1,nspinor 786 jspinor=3-ispinor 787 do ia=1,nincat 788 index_enl=atindx1(iatm+ia) 789 !$OMP DO 790 do jlmn=1,nlmn 791 j0lmn=jlmn*(jlmn-1)/2 792 jjlmn=j0lmn+jlmn 793 enl_(1:2)=enl(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor) 794 do mu=1,ndgxdtfac 795 if(cplex_dgxdt(mu)==2)then 796 cplex_ = 2 ; 797 gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,jlmn,ia,ispinor) 798 gxfj(1,mu) = zero ; gxfj(2,mu) = dgxdt(1,mu,jlmn,ia,jspinor) 799 else 800 cplex_ = cplex ; 801 gxfi(1:cplex) =dgxdt(1:cplex,mu,jlmn,ia,ispinor) 802 gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,jspinor) 803 end if 804 dgxdtfac(1,mu,jlmn,ia,jspinor)=dgxdtfac(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1) 805 dgxdtfac(2,mu,jlmn,ia,jspinor)=dgxdtfac(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1) 806 if (cplex_==2) then 807 dgxdtfac(2,mu,jlmn,ia,jspinor)=dgxdtfac(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2) 808 dgxdtfac(1,mu,jlmn,ia,jspinor)=dgxdtfac(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2) 809 end if 810 end do 811 do ilmn=1,jlmn-1 812 ijlmn=j0lmn+ilmn 813 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor) 814 do mu=1,ndgxdtfac 815 if(cplex_dgxdt(mu)==2)then 816 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,ispinor) 817 else 818 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 819 end if 820 dgxdtfac(1,mu,jlmn,ia,jspinor)=dgxdtfac(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1) 821 dgxdtfac(2,mu,jlmn,ia,jspinor)=dgxdtfac(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1) 822 #if !defined HAVE_OPENMP 823 dgxdtfac(1,mu,ilmn,ia,ispinor)=dgxdtfac(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu) 824 dgxdtfac(2,mu,ilmn,ia,ispinor)=dgxdtfac(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu) 825 #endif 826 if (cplex_==2) then 827 dgxdtfac(1,mu,jlmn,ia,jspinor)=dgxdtfac(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2) 828 dgxdtfac(2,mu,jlmn,ia,jspinor)=dgxdtfac(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2) 829 #if !defined HAVE_OPENMP 830 dgxdtfac(1,mu,ilmn,ia,ispinor)=dgxdtfac(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu) 831 dgxdtfac(2,mu,ilmn,ia,ispinor)=dgxdtfac(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu) 832 #endif 833 end if 834 end do !mu 835 end do !ilmn 836 #if defined HAVE_OPENMP 837 if(jlmn<nlmn) then 838 do ilmn=jlmn+1,nlmn 839 i0lmn=ilmn*(ilmn-1)/2 840 ijlmn=i0lmn+jlmn 841 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor) 842 do mu=1,ndgxdtfac 843 if(cplex_dgxdt(mu)==2)then 844 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,jspinor) 845 else 846 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,jspinor) 847 end if 848 dgxdtfac(1,mu,jlmn,ia,ispinor)=dgxdtfac(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1) 849 dgxdtfac(2,mu,jlmn,ia,ispinor)=dgxdtfac(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1) 850 if (cplex_==2) then 851 dgxdtfac(1,mu,jlmn,ia,ispinor)=dgxdtfac(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2) 852 dgxdtfac(2,mu,jlmn,ia,ispinor)=dgxdtfac(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2) 853 end if 854 end do !mu 855 end do !ilmn 856 end if 857 #endif 858 end do !jmln 859 !$OMP END DO 860 end do !ia 861 end do !ispinor 862 ABI_DEALLOCATE(gxfj) 863 !$OMP END PARALLEL 864 ! --- Parallelization over spinors --- 865 else if (nspinortot==2.and.nspinor/=nspinortot) then 866 ABI_CHECK(cplex_enl==2,"BUG: opernlc_ylm: invalid cplex_enl/=2!") 867 ABI_CHECK(cplex_fac==2,"BUG: opernlc_ylm: invalid cplex_fac/=2!") 868 ABI_ALLOCATE(dgxdtfac_,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinortot)) 869 !$OMP PARALLEL & 870 !$OMP PRIVATE(ia,index_enl), & 871 !$OMP PRIVATE(jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,mu,gxfi) 872 !$OMP WORKSHARE 873 dgxdtfac_(:,:,:,:,:)=zero 874 !$OMP END WORKSHARE 875 !$OMP SINGLE 876 ispinor_index=mpi_enreg%me_spinor+1 877 jspinor_index=3-ispinor_index 878 if (ispinor_index==1) then 879 ijspin=3;jispin=4 880 else 881 ijspin=4;jispin=3 882 end if 883 !$OMP END SINGLE 884 do ia=1,nincat 885 index_enl=atindx1(iatm+ia) 886 !$OMP DO 887 do jlmn=1,nlmn 888 j0lmn=jlmn*(jlmn-1)/2 889 do ilmn=1,nlmn 890 i0lmn=ilmn*(ilmn-1)/2 891 if (ilmn<=jlmn) then 892 ijlmn=j0lmn+ilmn 893 enl_(1)= enl(2*ijlmn-1,index_enl,ijspin) 894 enl_(2)=-enl(2*ijlmn ,index_enl,ijspin) 895 else 896 jilmn=i0lmn+jlmn 897 enl_(1:2)=enl(2*jilmn-1:2*jilmn,index_enl,jispin) 898 end if 899 do mu=1,ndgxdtfac 900 if(cplex_dgxdt(mu)==2)then 901 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = dgxdt(1,mu,ilmn,ia,1) 902 else 903 cplex_ = cplex ; gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,1) 904 end if 905 dgxdtfac_(1,mu,jlmn,ia,jspinor_index)=dgxdtfac_(1,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(1) 906 dgxdtfac_(2,mu,jlmn,ia,jspinor_index)=dgxdtfac_(2,mu,jlmn,ia,jspinor_index)+enl_(2)*gxfi(1) 907 if (cplex_==2) then 908 dgxdtfac_(1,mu,jlmn,ia,jspinor_index)=dgxdtfac_(1,mu,jlmn,ia,jspinor_index)-enl_(2)*gxfi(2) 909 dgxdtfac_(2,mu,jlmn,ia,jspinor_index)=dgxdtfac_(2,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(2) 910 end if 911 end do 912 end do !ilmn 913 end do !jlmn 914 !$OMP END DO 915 end do !iat 916 !$OMP SINGLE 917 call xmpi_sum(dgxdtfac_,mpi_enreg%comm_spinor,ierr) 918 !$OMP END SINGLE 919 !$OMP WORKSHARE 920 dgxdtfac(:,:,:,:,1)=dgxdtfac(:,:,:,:,1)+dgxdtfac_(:,:,:,:,ispinor_index) 921 !$OMP END WORKSHARE 922 !$OMP END PARALLEL 923 ABI_DEALLOCATE(dgxdtfac_) 924 end if !nspinortot 925 926 end if ! pawopt & optder 927 928 !Accumulate dgxdtfac related to overlap (Sij) (PAW) 929 !------------------------------------------------------------------- 930 if (optder>=1.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution 931 !$OMP PARALLEL & 932 !$OMP PRIVATE(ispinor,ia), & 933 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,sijr,mu,gxfj,ilmn,ijlmn,gxfi) 934 ABI_ALLOCATE(gxfj,(cplex,ndgxdtfac)) 935 !$OMP WORKSHARE 936 dgxdtfac_sij(1:cplex,1:ndgxdtfac,1:nlmn,1:nincat,1:nspinor)=zero 937 !$OMP END WORKSHARE 938 do ispinor=1,nspinor 939 do ia=1,nincat 940 !$OMP DO 941 do jlmn=1,nlmn 942 j0lmn=jlmn*(jlmn-1)/2 943 jjlmn=j0lmn+jlmn 944 sijr=sij(jjlmn) 945 do mu=1,ndgxdtfac 946 gxfj(1:cplex,mu)=dgxdt(1:cplex,mu,jlmn,ia,ispinor) 947 dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfj(1:cplex,mu) 948 end do 949 do ilmn=1,jlmn-1 950 ijlmn=j0lmn+ilmn 951 sijr=sij(ijlmn) 952 do mu=1,ndgxdtfac 953 gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 954 dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex) 955 #if !defined HAVE_OPENMP 956 dgxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)+sijr*gxfj(1:cplex,mu) 957 #endif 958 end do 959 end do 960 #if defined HAVE_OPENMP 961 if(jlmn<nlmn) then 962 do ilmn=jlmn+1,nlmn 963 i0lmn=ilmn*(ilmn-1)/2 964 ijlmn=i0lmn+jlmn 965 sijr=sij(ijlmn) 966 do mu=1,ndgxdtfac 967 gxfi(1:cplex)=dgxdt(1:cplex,mu,ilmn,ia,ispinor) 968 dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=dgxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex) 969 end do 970 end do 971 end if 972 #endif 973 end do 974 !$OMP END DO 975 end do 976 end do 977 ABI_DEALLOCATE(gxfj) 978 !$OMP END PARALLEL 979 end if 980 981 !Accumulate d2gxdtfac related to nonlocal operator (Norm-conserving) 982 !------------------------------------------------------------------- 983 if (optder==2.and.paw_opt==0) then ! Enl is E(Kleinman-Bylander) 984 !$OMP PARALLEL & 985 !$OMP PRIVATE(ispinor,ispinor_index,ia,ilmn,iln,enl_,mu) 986 do ispinor=1,nspinor 987 ispinor_index = ispinor + shift 988 do ia=1,nincat 989 !$OMP DO 990 do ilmn=1,nlmn 991 iln=indlmn(5,ilmn) 992 enl_(1)=enl(iln,itypat,ispinor_index) 993 do mu=1,nd2gxdtfac 994 d2gxdtfac(1:cplex,mu,ilmn,ia,ispinor)=enl_(1)*d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 995 end do 996 end do 997 !$OMP END DO 998 end do 999 end do 1000 !$OMP END PARALLEL 1001 end if 1002 1003 DBG_EXIT("COLL") 1004 1005 !Accumulate d2gxdtfac related to nonlocal operator (PAW) 1006 !------------------------------------------------------------------- 1007 if (optder==2.and.(paw_opt==1.or.paw_opt==2.or.paw_opt==4)) then ! Enl is psp strength Dij 1008 d2gxdtfac(1:cplex,1:nd2gxdtfac,1:nlmn,1:nincat,1:nspinor)=zero 1009 ! === Diagonal term(s) (up-up, down-down) 1010 ! 1-Enl is real 1011 if (cplex_enl==1) then 1012 !$OMP PARALLEL & 1013 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi) 1014 ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac)) 1015 do ispinor=1,nspinor 1016 ispinor_index=ispinor+shift 1017 do ia=1,nincat 1018 index_enl=atindx1(iatm+ia) 1019 !$OMP DO 1020 do jlmn=1,nlmn 1021 j0lmn=jlmn*(jlmn-1)/2 1022 jjlmn=j0lmn+jlmn 1023 enl_(1)=enl(jjlmn,index_enl,ispinor_index) 1024 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 1025 do mu=1,nd2gxdtfac 1026 gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor) 1027 d2gxdtfac(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu) 1028 end do 1029 do ilmn=1,jlmn-1 1030 ijlmn=j0lmn+ilmn 1031 enl_(1)=enl(ijlmn,index_enl,ispinor_index) 1032 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 1033 do mu=1,nd2gxdtfac 1034 gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1035 d2gxdtfac(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex) 1036 #if !defined HAVE_OPENMP 1037 d2gxdtfac(1:cplex,mu,ilmn,ia,ispinor)=d2gxdtfac(1:cplex,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1:cplex,mu) 1038 #endif 1039 end do 1040 end do 1041 #if defined HAVE_OPENMP 1042 if(jlmn<nlmn) then 1043 do ilmn=jlmn+1,nlmn 1044 i0lmn=ilmn*(ilmn-1)/2 1045 ijlmn=i0lmn+jlmn 1046 enl_(1)=enl(ijlmn,index_enl,ispinor_index) 1047 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 1048 do mu=1,nd2gxdtfac 1049 gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1050 d2gxdtfac(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac(1:cplex,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1:cplex) 1051 end do 1052 end do 1053 end if 1054 #endif 1055 end do 1056 !$OMP END DO 1057 end do 1058 end do 1059 ABI_DEALLOCATE(gxfj) 1060 !$OMP END PARALLEL 1061 ! 2-Enl is complex 1062 else 1063 ABI_CHECK(cplex_fac==cplex_enl,"BUG: invalid cplex_fac/=cplex_enl!") 1064 if (nspinortot==1) then !===== when nspinor=1, D_ij=D_ji 1065 !$OMP PARALLEL & 1066 !$OMP PRIVATE(ia,index_enl,jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi) 1067 ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac)) 1068 do ia=1,nincat 1069 index_enl=atindx1(iatm+ia) 1070 !$OMP DO 1071 do jlmn=1,nlmn 1072 j0lmn=jlmn*(jlmn-1)/2 1073 jjlmn=j0lmn+jlmn 1074 enl_(1:2)=enl(2*jjlmn-1:2*jjlmn,index_enl,1) 1075 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 1076 do mu=1,nd2gxdtfac 1077 if(cplex_d2gxdt(mu)==2)then 1078 cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,1) 1079 else 1080 cplex_ = cplex ; gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,1) 1081 end if 1082 d2gxdtfac(1,mu,jlmn,ia,1)=d2gxdtfac(1,mu,jlmn,ia,1)+enl_(1)*gxfj(1,mu) 1083 d2gxdtfac(2,mu,jlmn,ia,1)=d2gxdtfac(2,mu,jlmn,ia,1)+enl_(2)*gxfj(1,mu) 1084 if (cplex_==2) then 1085 d2gxdtfac(1,mu,jlmn,ia,1)=d2gxdtfac(1,mu,jlmn,ia,1)-enl_(2)*gxfj(2,mu) 1086 d2gxdtfac(2,mu,jlmn,ia,1)=d2gxdtfac(2,mu,jlmn,ia,1)+enl_(1)*gxfj(2,mu) 1087 end if 1088 end do 1089 do ilmn=1,jlmn-1 1090 ijlmn=j0lmn+ilmn 1091 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,1) 1092 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 1093 do mu=1,nd2gxdtfac 1094 if(cplex_d2gxdt(mu)==2)then 1095 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1) 1096 else 1097 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1) 1098 end if 1099 d2gxdtfac(1,mu,jlmn,ia,1)=d2gxdtfac(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1) 1100 d2gxdtfac(2,mu,jlmn,ia,1)=d2gxdtfac(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1) 1101 #if !defined HAVE_OPENMP 1102 d2gxdtfac(1,mu,ilmn,ia,1)=d2gxdtfac(1,mu,ilmn,ia,1)+enl_(1)*gxfj(1,mu) 1103 d2gxdtfac(2,mu,ilmn,ia,1)=d2gxdtfac(2,mu,ilmn,ia,1)+enl_(2)*gxfj(1,mu) 1104 #endif 1105 if (cplex_==2) then 1106 d2gxdtfac(1,mu,jlmn,ia,1)=d2gxdtfac(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2) 1107 d2gxdtfac(2,mu,jlmn,ia,1)=d2gxdtfac(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2) 1108 #if !defined HAVE_OPENMP 1109 d2gxdtfac(1,mu,ilmn,ia,1)=d2gxdtfac(1,mu,ilmn,ia,1)-enl_(2)*gxfj(2,mu) 1110 d2gxdtfac(2,mu,ilmn,ia,1)=d2gxdtfac(2,mu,ilmn,ia,1)+enl_(1)*gxfj(2,mu) 1111 #endif 1112 end if 1113 end do 1114 end do 1115 #if defined HAVE_OPENMP 1116 if(jlmn<nlmn) then 1117 do ilmn=jlmn+1,nlmn 1118 i0lmn=ilmn*(ilmn-1)/2 1119 ijlmn=i0lmn+jlmn 1120 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,1) 1121 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 1122 do mu=1,nd2gxdtfac 1123 if(cplex_d2gxdt(mu)==2)then 1124 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1) 1125 else 1126 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1) 1127 end if 1128 d2gxdtfac(1,mu,jlmn,ia,1)=d2gxdtfac(1,mu,jlmn,ia,1)+enl_(1)*gxfi(1) 1129 d2gxdtfac(2,mu,jlmn,ia,1)=d2gxdtfac(2,mu,jlmn,ia,1)+enl_(2)*gxfi(1) 1130 if (cplex_==2) then 1131 d2gxdtfac(1,mu,jlmn,ia,1)=d2gxdtfac(1,mu,jlmn,ia,1)-enl_(2)*gxfi(2) 1132 d2gxdtfac(2,mu,jlmn,ia,1)=d2gxdtfac(2,mu,jlmn,ia,1)+enl_(1)*gxfi(2) 1133 end if 1134 end do 1135 end do 1136 end if 1137 #endif 1138 end do 1139 !$OMP END DO 1140 end do 1141 ABI_DEALLOCATE(gxfj) 1142 !$OMP END PARALLEL 1143 else !===== when nspinor=2, D_ij=D_ji^* 1144 !$OMP PARALLEL & 1145 !$OMP PRIVATE(ispinor,ispinor_index,ia,index_enl), & 1146 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfj,ilmn,ijlmn,gxfi) 1147 ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac)) 1148 do ispinor=1,nspinor 1149 ispinor_index = ispinor + shift 1150 do ia=1,nincat 1151 index_enl=atindx1(iatm+ia) 1152 !$OMP DO 1153 do jlmn=1,nlmn 1154 j0lmn=jlmn*(jlmn-1)/2 1155 jjlmn=j0lmn+jlmn 1156 enl_(1)=enl(2*jjlmn-1,index_enl,ispinor_index) ! enl_ii is real for up-up and dn-dn 1157 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(jjlmn) 1158 do mu=1,nd2gxdtfac 1159 if(cplex_d2gxdt(mu)==2)then 1160 cplex_ = 2 ; gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,ispinor) 1161 else 1162 cplex_ = cplex ; gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor) 1163 end if 1164 d2gxdtfac(1:cplex_,mu,jlmn,ia,ispinor)=d2gxdtfac(1:cplex_,mu,jlmn,ia,ispinor)+enl_(1)*gxfj(1:cplex_,mu) 1165 end do 1166 do ilmn=1,jlmn-1 1167 ijlmn=j0lmn+ilmn 1168 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index) 1169 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 1170 do mu=1,nd2gxdtfac 1171 if(cplex_d2gxdt(mu)==2)then 1172 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor) 1173 else 1174 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1175 end if 1176 d2gxdtfac(1,mu,jlmn,ia,ispinor)=d2gxdtfac(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1) 1177 d2gxdtfac(2,mu,jlmn,ia,ispinor)=d2gxdtfac(2,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(1) 1178 #if !defined HAVE_OPENMP 1179 d2gxdtfac(1,mu,ilmn,ia,ispinor)=d2gxdtfac(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu) 1180 d2gxdtfac(2,mu,ilmn,ia,ispinor)=d2gxdtfac(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu) 1181 #endif 1182 if (cplex_==2) then 1183 d2gxdtfac(1,mu,jlmn,ia,ispinor)=d2gxdtfac(1,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(2) 1184 d2gxdtfac(2,mu,jlmn,ia,ispinor)=d2gxdtfac(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2) 1185 #if !defined HAVE_OPENMP 1186 d2gxdtfac(1,mu,ilmn,ia,ispinor)=d2gxdtfac(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu) 1187 d2gxdtfac(2,mu,ilmn,ia,ispinor)=d2gxdtfac(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu) 1188 #endif 1189 end if 1190 end do 1191 end do 1192 #if defined HAVE_OPENMP 1193 if(jlmn<nlmn) then 1194 do ilmn=jlmn+1,nlmn 1195 i0lmn=ilmn*(ilmn-1)/2 1196 ijlmn=i0lmn+jlmn 1197 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,ispinor_index) 1198 if (paw_opt==2) enl_(1)=enl_(1)-lambda*sij(ijlmn) 1199 do mu=1,nd2gxdtfac 1200 if(cplex_d2gxdt(mu)==2)then 1201 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor) 1202 else 1203 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1204 end if 1205 d2gxdtfac(1,mu,jlmn,ia,ispinor)=d2gxdtfac(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1) 1206 d2gxdtfac(2,mu,jlmn,ia,ispinor)=d2gxdtfac(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1) 1207 if (cplex_==2) then 1208 d2gxdtfac(1,mu,jlmn,ia,ispinor)=d2gxdtfac(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2) 1209 d2gxdtfac(2,mu,jlmn,ia,ispinor)=d2gxdtfac(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2) 1210 end if 1211 end do 1212 end do 1213 end if 1214 #endif 1215 end do 1216 !$OMP END DO 1217 end do 1218 end do 1219 ABI_DEALLOCATE(gxfj) 1220 !$OMP END PARALLEL 1221 end if !nspinortot 1222 end if !complex 1223 ! === Off-diagonal term(s) (up-down, down-up) 1224 ! --- No parallelization over spinors --- 1225 if (nspinortot==2.and.nspinor==nspinortot) then 1226 ABI_CHECK(cplex_enl==2,"BUG: invalid cplex_enl/=2!") 1227 ABI_CHECK(cplex_fac==2,"BUG: invalid cplex_fac/=2!") 1228 !$OMP PARALLEL & 1229 !$OMP PRIVATE(ispinor,jspinor,ia,index_enl), & 1230 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,enl_,mu,gxfi,gxfj,ilmn,ijlmn) 1231 ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac)) 1232 do ispinor=1,nspinor 1233 jspinor=3-ispinor 1234 do ia=1,nincat 1235 index_enl=atindx1(iatm+ia) 1236 !$OMP DO 1237 do jlmn=1,nlmn 1238 j0lmn=jlmn*(jlmn-1)/2 1239 jjlmn=j0lmn+jlmn 1240 enl_(1:2)=enl(2*jjlmn-1:2*jjlmn,index_enl,2+ispinor) 1241 do mu=1,nd2gxdtfac 1242 if(cplex_d2gxdt(mu)==2)then 1243 cplex_ = 2 ; 1244 gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,jlmn,ia,ispinor) 1245 gxfj(1,mu) = zero ; gxfj(2,mu) = d2gxdt(1,mu,jlmn,ia,jspinor) 1246 else 1247 cplex_ = cplex ; 1248 gxfi(1:cplex) =d2gxdt(1:cplex,mu,jlmn,ia,ispinor) 1249 gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,jspinor) 1250 end if 1251 d2gxdtfac(1,mu,jlmn,ia,jspinor)=d2gxdtfac(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1) 1252 d2gxdtfac(2,mu,jlmn,ia,jspinor)=d2gxdtfac(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1) 1253 if (cplex_==2) then 1254 d2gxdtfac(2,mu,jlmn,ia,jspinor)=d2gxdtfac(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2) 1255 d2gxdtfac(1,mu,jlmn,ia,jspinor)=d2gxdtfac(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2) 1256 end if 1257 end do 1258 do ilmn=1,jlmn-1 1259 ijlmn=j0lmn+ilmn 1260 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor) 1261 do mu=1,nd2gxdtfac 1262 if(cplex_d2gxdt(mu)==2)then 1263 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,ispinor) 1264 else 1265 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1266 end if 1267 d2gxdtfac(1,mu,jlmn,ia,jspinor)=d2gxdtfac(1,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(1) 1268 d2gxdtfac(2,mu,jlmn,ia,jspinor)=d2gxdtfac(2,mu,jlmn,ia,jspinor)-enl_(2)*gxfi(1) 1269 #if !defined HAVE_OPENMP 1270 d2gxdtfac(1,mu,ilmn,ia,ispinor)=d2gxdtfac(1,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(1,mu) 1271 d2gxdtfac(2,mu,ilmn,ia,ispinor)=d2gxdtfac(2,mu,ilmn,ia,ispinor)+enl_(2)*gxfj(1,mu) 1272 #endif 1273 if (cplex_==2) then 1274 d2gxdtfac(1,mu,jlmn,ia,jspinor)=d2gxdtfac(1,mu,jlmn,ia,jspinor)+enl_(2)*gxfi(2) 1275 d2gxdtfac(2,mu,jlmn,ia,jspinor)=d2gxdtfac(2,mu,jlmn,ia,jspinor)+enl_(1)*gxfi(2) 1276 #if !defined HAVE_OPENMP 1277 d2gxdtfac(1,mu,ilmn,ia,ispinor)=d2gxdtfac(1,mu,ilmn,ia,ispinor)-enl_(2)*gxfj(2,mu) 1278 d2gxdtfac(2,mu,ilmn,ia,ispinor)=d2gxdtfac(2,mu,ilmn,ia,ispinor)+enl_(1)*gxfj(2,mu) 1279 #endif 1280 end if 1281 end do !mu 1282 end do !ilmn 1283 #if defined HAVE_OPENMP 1284 if(jlmn<nlmn) then 1285 do ilmn=jlmn+1,nlmn 1286 i0lmn=ilmn*(ilmn-1)/2 1287 ijlmn=i0lmn+jlmn 1288 enl_(1:2)=enl(2*ijlmn-1:2*ijlmn,index_enl,2+ispinor) 1289 do mu=1,nd2gxdtfac 1290 if(cplex_d2gxdt(mu)==2)then 1291 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,jspinor) 1292 else 1293 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,jspinor) 1294 end if 1295 d2gxdtfac(1,mu,jlmn,ia,ispinor)=d2gxdtfac(1,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(1) 1296 d2gxdtfac(2,mu,jlmn,ia,ispinor)=d2gxdtfac(2,mu,jlmn,ia,ispinor)+enl_(2)*gxfi(1) 1297 if (cplex_==2) then 1298 d2gxdtfac(1,mu,jlmn,ia,ispinor)=d2gxdtfac(1,mu,jlmn,ia,ispinor)-enl_(2)*gxfi(2) 1299 d2gxdtfac(2,mu,jlmn,ia,ispinor)=d2gxdtfac(2,mu,jlmn,ia,ispinor)+enl_(1)*gxfi(2) 1300 end if 1301 end do !mu 1302 end do !ilmn 1303 end if 1304 #endif 1305 end do !jmln 1306 !$OMP END DO 1307 end do !ia 1308 end do !ispinor 1309 ABI_DEALLOCATE(gxfj) 1310 !$OMP END PARALLEL 1311 ! --- Parallelization over spinors --- 1312 else if (nspinortot==2.and.nspinor/=nspinortot) then 1313 ABI_CHECK(cplex_enl==2,"BUG: opernlc_ylm: invalid cplex_enl/=2!") 1314 ABI_CHECK(cplex_fac==2,"BUG: opernlc_ylm: invalid cplex_fac/=2!") 1315 ABI_ALLOCATE(d2gxdtfac_,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinortot)) 1316 !$OMP PARALLEL & 1317 !$OMP PRIVATE(ia,index_enl), & 1318 !$OMP PRIVATE(jlmn,j0lmn,ilmn,i0lmn,ijlmn,enl_,jilmn,mu,gxfi) 1319 !$OMP WORKSHARE 1320 d2gxdtfac_(:,:,:,:,:)=zero 1321 !$OMP END WORKSHARE 1322 !$OMP SINGLE 1323 ispinor_index=mpi_enreg%me_spinor+1 1324 jspinor_index=3-ispinor_index 1325 if (ispinor_index==1) then 1326 ijspin=3;jispin=4 1327 else 1328 ijspin=4;jispin=3 1329 end if 1330 !$OMP END SINGLE 1331 do ia=1,nincat 1332 index_enl=atindx1(iatm+ia) 1333 !$OMP DO 1334 do jlmn=1,nlmn 1335 j0lmn=jlmn*(jlmn-1)/2 1336 do ilmn=1,nlmn 1337 i0lmn=ilmn*(ilmn-1)/2 1338 if (ilmn<=jlmn) then 1339 ijlmn=j0lmn+ilmn 1340 enl_(1)= enl(2*ijlmn-1,index_enl,ijspin) 1341 enl_(2)=-enl(2*ijlmn ,index_enl,ijspin) 1342 else 1343 jilmn=i0lmn+jlmn 1344 enl_(1:2)=enl(2*jilmn-1:2*jilmn,index_enl,jispin) 1345 end if 1346 do mu=1,nd2gxdtfac 1347 if(cplex_d2gxdt(mu)==2)then 1348 cplex_ = 2 ; gxfi(1) = zero ; gxfi(2) = d2gxdt(1,mu,ilmn,ia,1) 1349 else 1350 cplex_ = cplex ; gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,1) 1351 end if 1352 d2gxdtfac_(1,mu,jlmn,ia,jspinor_index)=d2gxdtfac_(1,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(1) 1353 d2gxdtfac_(2,mu,jlmn,ia,jspinor_index)=d2gxdtfac_(2,mu,jlmn,ia,jspinor_index)+enl_(2)*gxfi(1) 1354 if (cplex_==2) then 1355 d2gxdtfac_(1,mu,jlmn,ia,jspinor_index)=d2gxdtfac_(1,mu,jlmn,ia,jspinor_index)-enl_(2)*gxfi(2) 1356 d2gxdtfac_(2,mu,jlmn,ia,jspinor_index)=d2gxdtfac_(2,mu,jlmn,ia,jspinor_index)+enl_(1)*gxfi(2) 1357 end if 1358 end do 1359 end do !ilmn 1360 end do !jlmn 1361 !$OMP END DO 1362 end do !iat 1363 !$OMP SINGLE 1364 call xmpi_sum(d2gxdtfac_,mpi_enreg%comm_spinor,ierr) 1365 !$OMP END SINGLE 1366 !$OMP WORKSHARE 1367 d2gxdtfac(:,:,:,:,1)=d2gxdtfac(:,:,:,:,1)+d2gxdtfac_(:,:,:,:,ispinor_index) 1368 !$OMP END WORKSHARE 1369 !$OMP END PARALLEL 1370 ABI_DEALLOCATE(d2gxdtfac_) 1371 end if !nspinortot 1372 1373 end if ! pawopt & optder 1374 1375 !Accumulate d2gxdtfac related to overlap (Sij) (PAW) 1376 !------------------------------------------------------------------- 1377 if (optder==2.and.(paw_opt==3.or.paw_opt==4)) then ! Use Sij, overlap contribution 1378 !$OMP PARALLEL & 1379 !$OMP PRIVATE(ispinor,ia), & 1380 !$OMP PRIVATE(jlmn,j0lmn,jjlmn,sijr,mu,gxfj,ilmn,ijlmn,gxfi) 1381 ABI_ALLOCATE(gxfj,(cplex,nd2gxdtfac)) 1382 !$OMP WORKSHARE 1383 d2gxdtfac_sij(1:cplex,1:nd2gxdtfac,1:nlmn,1:nincat,1:nspinor)=zero 1384 !$OMP END WORKSHARE 1385 do ispinor=1,nspinor 1386 do ia=1,nincat 1387 !$OMP DO 1388 do jlmn=1,nlmn 1389 j0lmn=jlmn*(jlmn-1)/2 1390 jjlmn=j0lmn+jlmn 1391 sijr=sij(jjlmn) 1392 do mu=1,nd2gxdtfac 1393 gxfj(1:cplex,mu)=d2gxdt(1:cplex,mu,jlmn,ia,ispinor) 1394 d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfj(1:cplex,mu) 1395 end do 1396 do ilmn=1,jlmn-1 1397 ijlmn=j0lmn+ilmn 1398 sijr=sij(ijlmn) 1399 do mu=1,nd2gxdtfac 1400 gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1401 d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex) 1402 #if !defined HAVE_OPENMP 1403 d2gxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)=d2gxdtfac_sij(1:cplex,mu,ilmn,ia,ispinor)+sijr*gxfj(1:cplex,mu) 1404 #endif 1405 end do 1406 end do 1407 #if defined HAVE_OPENMP 1408 if(jlmn<nlmn) then 1409 do ilmn=jlmn+1,nlmn 1410 i0lmn=ilmn*(ilmn-1)/2 1411 ijlmn=i0lmn+jlmn 1412 sijr=sij(ijlmn) 1413 do mu=1,nd2gxdtfac 1414 gxfi(1:cplex)=d2gxdt(1:cplex,mu,ilmn,ia,ispinor) 1415 d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)=d2gxdtfac_sij(1:cplex,mu,jlmn,ia,ispinor)+sijr*gxfi(1:cplex) 1416 end do 1417 end do 1418 end if 1419 #endif 1420 end do 1421 !$OMP END DO 1422 end do 1423 end do 1424 ABI_DEALLOCATE(gxfj) 1425 !$OMP END PARALLEL 1426 end if 1427 1428 end subroutine opernlc_ylm