TABLE OF CONTENTS
ABINIT/opernld_ylm [ Functions ]
NAME
opernld_ylm
FUNCTION
* Operate with the non-local part of the hamiltonian, in order to get contributions to energy/forces/stress/dyn.matrix/elst tens. from projected scalars * Operate with the non-local projectors and the overlap matrix Sij in order to get contributions to <c|S|c> from 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
choice=chooses possible output
cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
2 if <p_lmn|c> scalars are complex
cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)=2nd gradients of projected scalars
dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)=gradients of projected scalars
dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
related to Vnl (NL operator)
dgxdtfac_sij(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
related to Sij (overlap)
gx(cplex,nlmn,nincat,nspinor)= projected scalars
gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
ia3=gives the absolute number of the first atom in the subset presently treated
natom=number of atoms in cell
nd2gxdt=second dimension of d2gxdt
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
nnlout=dimension of enlout
nspinor=number of spinorial components of the wavefunctions (on current proc)
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)
OUTPUT
(see side effects)
SIDE EFFECTS
--If (paw_opt==0, 1 or 2)
enlout(nnlout)= contribution to the non-local part of the following properties:
if choice=1 : enlout(1) -> the energy
if choice=2 : enlout(3*natom) -> 1st deriv. of energy wrt atm. pos (forces)
if choice=3 : enlout(6) -> 1st deriv. of energy wrt strain (stresses)
if choice=4 : enlout(6*natom) -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
if choice=23: enlout(6+3*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
1st deriv. of energy wrt strain (stresses)
if choice=24: enlout(9*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
if choice=5 : enlout(3) -> 1st deriv. of energy wrt k
if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
if choice=6 : enlout(36+18*natom) -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
2nd deriv. of energy wrt to atm. pos and strain (internal strain)
if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
if choice=81: enlout(18) -> 2nd deriv. of energy wrt k and right k
--If (paw_opt==3)
if choice=1 : enlout(1) -> contribution to <c|S|c> (note: not including <c|c>)
if choice=2 : enlout(3*natom) -> contribution to <c|dS/d_atm.pos|c>
if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
if choice=81: enlout(18) -> 2nd deriv. of energy wrt k and right k
--If (paw_opt==4)
not available
NOTES
Operate for one type of atom, and within this given type of atom, for a subset of at most nincat atoms.
PARENTS
nonlop_ylm
CHILDREN
SOURCE
94 #if defined HAVE_CONFIG_H 95 #include "config.h" 96 #endif 97 98 #include "abi_common.h" 99 100 101 subroutine opernld_ylm(choice,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,& 102 & enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,nd2gxdt,ndgxdt,& 103 & ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt,strnlk) 104 105 106 use defs_basis 107 use m_profiling_abi 108 use m_errors 109 110 !This section has been created automatically by the script Abilint (TD). 111 !Do not modify the following lines by hand. 112 #undef ABI_FUNC 113 #define ABI_FUNC 'opernld_ylm' 114 !End of the abilint section 115 116 implicit none 117 118 !Arguments ------------------------------------ 119 !scalars 120 integer,intent(in) :: choice,cplex,cplex_fac,ia3,natom,nd2gxdt,ndgxdt 121 integer,intent(in) :: ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt 122 real(dp),intent(inout) :: enlk 123 !arrays 124 real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor) 125 real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor) 126 real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor) 127 real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3)) 128 real(dp),intent(in) :: gx(cplex,nlmn,nincat,nspinor),gxfac(cplex_fac,nlmn,nincat,nspinor) 129 real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3)) 130 real(dp),intent(inout) :: ddkk(6),enlout(nnlout),fnlk(3*natom),strnlk(6) 131 132 !Local variables------------------------------- 133 !scalars 134 integer :: ia,iashift,ilmn,iplex,ishift,ispinor,mu,mua,mua1,mua2,mub,mushift,mut,muu,nu,nushift 135 real(dp) :: dummy 136 !arrays 137 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 138 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 139 integer,parameter :: twist_dir(6)=(/2,3,3,1,1,2/) 140 real(dp) :: d2gx(cplex),enlj(6*cplex),gxfacj(cplex) 141 real(dp),allocatable :: enljj(:) 142 complex(dpc),allocatable :: cft(:,:), cfu(:,:) 143 144 ! ************************************************************************* 145 146 ABI_CHECK(cplex_fac>=cplex,'BUG: invalid cplex_fac<cplex!') 147 148 if (paw_opt==0.or.paw_opt==1.or.paw_opt==2) then 149 150 ! ============== Accumulate the non-local energy =============== 151 if (choice==1) then 152 do ispinor=1,nspinor 153 do ia=1,nincat 154 do ilmn=1,nlmn 155 do iplex=1,cplex 156 enlout(1)=enlout(1)+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor) 157 end do 158 end do 159 end do 160 end do 161 end if 162 163 ! ============ Accumulate the forces contributions ============= 164 if (choice==2.or.choice==23.or.choice==24) then 165 ishift=0;if (choice==23) ishift=6 166 do ispinor=1,nspinor 167 do ia=1,nincat 168 enlj(1:3)=zero 169 iashift=3*(ia+ia3-2)+ishift 170 do ilmn=1,nlmn 171 do mu=1,3 172 dummy = zero ! Dummy needed here to get the correct forces with intel -O3 173 do iplex=1,cplex 174 !enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor) 175 dummy=dummy+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor) 176 end do 177 enlj(mu)=enlj(mu)+dummy 178 end do 179 end do 180 enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3) 181 end do 182 end do 183 end if 184 185 ! ======== Accumulate the stress tensor contributions ========== 186 if (choice==3.or.choice==23) then 187 enlj(1:6)=zero 188 do ispinor=1,nspinor 189 do ia=1,nincat 190 do ilmn=1,nlmn 191 gxfacj(1:cplex)=gxfac(1:cplex,ilmn,ia,ispinor) 192 do iplex=1,cplex 193 enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor) 194 end do 195 do mu=1,6 196 do iplex=1,cplex 197 enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor) 198 end do 199 end do 200 end do 201 end do 202 end do 203 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 204 end if 205 206 ! ====== Accumulate the dynamical matrix contributions ========= 207 if (choice==4.or.choice==24) then 208 ishift=0;if (choice==24) ishift=3*natom 209 do ispinor=1,nspinor 210 do ia=1,nincat 211 enlj(1:6)=zero 212 iashift=6*(ia+ia3-2)+ishift 213 do ilmn=1,nlmn 214 do mu=1,6 215 mua=alpha(mu);mub=beta(mu) 216 do iplex=1,cplex 217 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)& 218 & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,mua,ilmn,ia,ispinor) 219 end do 220 end do 221 end do 222 enlout(iashift+1:iashift+6)=enlout(iashift+1:iashift+6)+two*enlj(1:6) 223 end do 224 end do 225 end if 226 227 ! ======== Accumulate the contributions of derivatives of E wrt to k ========== 228 if (choice==5) then 229 enlj(1:3)=zero 230 do ispinor=1,nspinor 231 do ia=1,nincat 232 if(cplex==2)then 233 do ilmn=1,nlmn 234 do mu=1,3 235 do iplex=1,cplex 236 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor) 237 end do 238 end do 239 end do 240 ! If cplex=1, dgxdt is pure imaginary; thus there is no contribution 241 else if (cplex_fac==2) then 242 do ilmn=1,nlmn 243 do mu=1,3 244 enlj(mu)=enlj(mu)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 245 end do 246 end do 247 end if 248 end do 249 end do 250 enlout(1:3)=enlout(1:3)+two*enlj(1:3) 251 end if 252 253 ! ======== Accumulate the contributions of partial derivatives of E wrt to k ========== 254 ! Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k 255 if (choice==51.or.choice==52) then 256 enlj(1:6)=zero 257 do ispinor=1,nspinor 258 do ia=1,nincat 259 if(cplex==2)then 260 do ilmn=1,nlmn 261 do mu=1,3 262 enlj(2*mu-1)=enlj(2*mu-1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) & 263 & +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) 264 enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) & 265 & -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 266 end do 267 end do 268 else if (cplex_fac==2) then 269 do ilmn=1,nlmn 270 do mu=1,3 271 enlj(2*mu-1)=enlj(2*mu-1)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 272 enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 273 end do 274 end do 275 else if (cplex_fac==1) then 276 do ilmn=1,nlmn 277 do mu=1,3 278 enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 279 end do 280 end do 281 end if 282 end do 283 end do 284 if (choice==52) then 285 enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6) 286 end if 287 enlout(1:6)=enlout(1:6)+enlj(1:6) 288 end if 289 290 ! ======== Accumulate the contributions of twist derivatives of E wrt to k ========== 291 if (choice==53) then 292 enlj(1:3)=zero 293 ABI_ALLOCATE(cft,(3,nlmn)) 294 ABI_ALLOCATE(cfu,(3,nlmn)) 295 ! If cplex=1, dgxdt is pure imaginary; 296 ! If cplex_fac=1, dgxdtfac is pure imaginary; 297 do ispinor=1,nspinor 298 do ia=1,nincat 299 if(cplex==2)then 300 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 301 else 302 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 303 end if 304 if(cplex_fac==2)then 305 cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor)) 306 else 307 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor)) 308 end if 309 do ilmn=1,nlmn 310 do mu=1,3 311 mut = twist_dir(2*mu-1) 312 muu = twist_dir(2*mu) 313 enlj(mu) = enlj(mu) + aimag(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 314 enlj(mu) = enlj(mu) - aimag(conjg(cft(muu,ilmn))*cfu(mut,ilmn)) 315 end do 316 end do 317 end do 318 end do 319 enlout(1:3)=enlout(1:3)+enlj(1:3) 320 ABI_DEALLOCATE(cft) 321 ABI_DEALLOCATE(cfu) 322 end if 323 324 ! ====== Accumulate the effective charges contributions ========= 325 if (choice==54) then 326 ABI_ALLOCATE(enljj,(18)) 327 do ispinor=1,nspinor 328 do ia=1,nincat 329 enljj(1:18)=zero 330 iashift=18*(ia+ia3-2) 331 ! If cplex=1, dgxdt is real for atm. pos, pure imaginary for k; 332 ! If cplex_fac=1, dgxdtfac is pure imaginary for k; 333 if(cplex==2.and.cplex_fac==2) then 334 do ilmn=1,nlmn 335 mu=1;nu=1 336 do mua=1,3 ! atm. pos 337 do mub=1,3 ! k 338 enljj(nu)=enljj(nu) & 339 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 340 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) & 341 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) & 342 & +gxfac(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) 343 enljj(nu+1)=enljj(nu+1) & 344 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) & 345 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 346 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) & 347 & -gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 348 mu=mu+1;nu=nu+2 349 end do 350 end do 351 end do 352 else if(cplex==1.and.cplex_fac==2)then 353 do ilmn=1,nlmn 354 mu=1;nu=1 355 do mua=1,3 ! atm. pos 356 do mub=1,3 ! k 357 enljj(nu)=enljj(nu) & 358 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 359 & +gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 360 enljj(nu+1)=enljj(nu+1) & 361 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) & 362 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 363 mu=mu+1;nu=nu+2 364 end do 365 end do 366 end do 367 else if(cplex==1.and.cplex_fac==1)then 368 do ilmn=1,nlmn 369 mu=1;nu=1 370 do mua=1,3 ! atm. pos 371 do mub=1,3 ! k 372 enljj(nu+1)=enljj(nu+1) & 373 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 374 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 375 mu=mu+1;nu=nu+2 376 end do 377 end do 378 end do 379 end if 380 enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18) 381 end do 382 end do 383 ABI_DEALLOCATE(enljj) 384 end if 385 386 ! ====== Accumulate the piezoelectric tensor contributions ========= 387 if (choice==55) then 388 ABI_ALLOCATE(enljj,(36)) 389 do ispinor=1,nspinor 390 do ia=1,nincat 391 enljj(1:36)=zero;enlj(:)=zero 392 ! If cplex=1, dgxdt is real for strain, pure imaginary for k; 393 ! If cplex_fac=1, dgxdtfac is pure imaginary for k; 394 if(cplex==2.and.cplex_fac==2) then 395 do ilmn=1,nlmn 396 ! First compute 2nd-derivative contribution 397 mu=1 398 do mua=1,6 ! strain (lambda,nu) 399 mua1=alpha(mua) ! (nu) 400 mua2=beta(mua) ! (lambda) 401 do mub=1,3 ! k (mu) 402 muu=3*(gamma(mua1,mub)-1)+mua2 403 mut=3*(gamma(mua2,mub)-1)+mua1 404 d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) & 405 & +d2gxdt(1:cplex,mut,ilmn,ia,ispinor)) 406 enljj(mu)=enljj(mu) & 407 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 408 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) & 409 & +gxfac(1,ilmn,ia,ispinor)*d2gx(1)+gxfac(2,ilmn,ia,ispinor)*d2gx(2) 410 enljj(mu+1)=enljj(mu+1) & 411 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) & 412 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 413 & +gxfac(1,ilmn,ia,ispinor)*d2gx(2)-gxfac(2,ilmn,ia,ispinor)*d2gx(1) 414 mu=mu+2 415 end do 416 end do 417 ! Then store 1st-derivative contribution 418 mu=1 419 do nu=1,3 420 enlj(mu )=enlj(mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) & 421 & +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) 422 enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) & 423 & -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 424 mu=mu+2 425 end do 426 end do 427 else if(cplex==1.and.cplex_fac==2)then 428 do ilmn=1,nlmn 429 ! First compute 2nd-derivative contribution 430 mu=1 431 do mua=1,6 ! strain (lambda,nu) 432 mua1=alpha(mua) ! (nu) 433 mua2=beta(mua) ! (lambda) 434 do mub=1,3 ! k (mu) 435 muu=3*(gamma(mua1,mub)-1)+mua2 436 mut=3*(gamma(mua2,mub)-1)+mua1 437 d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor)) 438 enljj(mu)=enljj(mu) & 439 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 440 & +gxfac(2,ilmn,ia,ispinor)*d2gx(1) 441 enljj(mu+1)=enljj(mu+1) & 442 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) & 443 & +gxfac(1,ilmn,ia,ispinor)*d2gx(1) 444 mu=mu+2 445 end do 446 end do 447 ! Then store 1st-derivative contribution 448 mu=1 449 do nu=1,3 450 enlj(mu )=enlj(mu )+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 451 enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 452 mu=mu+2 453 end do 454 end do 455 else if(cplex==1.and.cplex_fac==1)then 456 do ilmn=1,nlmn 457 mu=1 458 do mua=1,6 ! strain (lambda,nu) 459 mua1=alpha(mua) ! (nu) 460 mua2=beta(mua) ! (lambda) 461 do mub=1,3 ! k (mu) 462 muu=3*(gamma(mua1,mub)-1)+mua2 463 mut=3*(gamma(mua2,mub)-1)+mua1 464 d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor)) 465 enljj(mu+1)=enljj(mu+1) & 466 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 467 & +gxfac(1,ilmn,ia,ispinor)*d2gx(1) 468 mu=mu+2 469 end do 470 end do 471 ! Then store 1st-derivative contribution 472 mu=1 473 do nu=1,3 474 enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 475 mu=mu+2 476 end do 477 end do 478 end if 479 enlout(1:36)=enlout(1:36)+enljj(1:36) 480 ddkk(1:6)=ddkk(1:6)+enlj(1:6) 481 end do 482 end do 483 ABI_DEALLOCATE(enljj) 484 end if 485 486 ! ======= Accumulate the elastic tensor contributions ========== 487 if (choice==6) then 488 do ispinor=1,nspinor 489 do ia=1,nincat 490 iashift=3*(ia+ia3-2) 491 do ilmn=1,nlmn 492 do iplex=1,cplex 493 enlk=enlk+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor) 494 end do 495 enlj(1:3)=zero 496 do mu=1,3 497 do iplex=1,cplex 498 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,6+mu,ilmn,ia,ispinor) 499 end do 500 end do 501 fnlk(iashift+1:iashift+3)=fnlk(iashift+1:iashift+3)+two*enlj(1:3) 502 enlj(1:6)=zero 503 do mu=1,6 504 do iplex=1,cplex 505 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor) 506 end do 507 end do 508 strnlk(1:6)=strnlk(1:6)+two*enlj(1:6) 509 do mub=1,6 510 mushift=6*(mub-1);nushift=(3*natom+6)*(mub-1) 511 do mua=1,6 512 mu=mushift+mua;nu=nushift+mua 513 do iplex=1,cplex 514 enlout(nu)=enlout(nu)+two* & 515 & (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)& 516 & +dgxdtfac(iplex,mua,ilmn,ia,ispinor)*dgxdt(iplex,mub,ilmn,ia,ispinor)) 517 end do 518 end do 519 mushift=36+3*(mub-1);nushift=6+iashift+(3*natom+6)*(mub-1) 520 do mua=1,3 521 mu=mushift+mua;nu=nushift+mua 522 do iplex=1,cplex 523 enlout(nu)=enlout(nu)+two* & 524 & (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)& 525 & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,6+mua,ilmn,ia,ispinor)) 526 end do 527 end do 528 end do 529 end do 530 end do 531 end do 532 end if 533 534 ! ======== Accumulate the contributions of 2nd-derivatives of E wrt to k ========== 535 if (choice==8) then 536 ABI_ALLOCATE(cft,(3,nlmn)) 537 ABI_ALLOCATE(cfu,(3,nlmn)) 538 do ispinor=1,nspinor 539 do ia=1,nincat 540 enlj(1:6)=zero 541 do ilmn=1,nlmn 542 do mu=1,6 543 do iplex=1,cplex 544 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor) 545 end do 546 end do 547 end do 548 ! If cplex=1, dgxdt is pure imaginary; 549 ! If cplex_fac=1, dgxdtfac is pure imaginary; 550 if(cplex==2)then 551 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 552 else 553 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 554 end if 555 if(cplex_fac==2)then 556 cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor)) 557 else 558 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor)) 559 end if 560 do ilmn=1,nlmn 561 do mu=1,6 562 mua=alpha(mu);mub=beta(mu) 563 enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn)) 564 end do 565 end do 566 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 567 end do 568 end do 569 ABI_DEALLOCATE(cft) 570 ABI_DEALLOCATE(cfu) 571 end if 572 573 ! ======== Accumulate the contributions of partial 2nd-derivatives of E wrt to k ========== 574 ! Full derivative wrt to k1, right derivative wrt to k2 575 if (choice==81) then 576 ABI_ALLOCATE(cft,(3,nlmn)) 577 ABI_ALLOCATE(cfu,(6,nlmn)) 578 ABI_ALLOCATE(enljj,(18)) 579 do ispinor=1,nspinor 580 do ia=1,nincat 581 enljj(1:18)=zero 582 if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real 583 cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),gxfac(2,1:nlmn,ia,ispinor)) 584 else 585 cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),zero) 586 end if 587 if(cplex==2)then !If cplex=1, d2gxdt is pure real 588 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor)) 589 else 590 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero) 591 end if 592 do ilmn=1,nlmn 593 do mu=1,3 594 do nu=1,3 595 muu=3*(mu-1)+nu ; mut=gamma(mu,nu) 596 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 597 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 598 end do 599 end do 600 end do 601 if(cplex==2)then !If cplex=1, dgxdt is pure imaginary 602 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 603 else 604 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 605 end if 606 if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary 607 cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor)) 608 else 609 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor)) 610 end if 611 do ilmn=1,nlmn 612 do mu=1,3 613 do nu=1,3 614 muu=3*(mu-1)+nu 615 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 616 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 617 end do 618 end do 619 end do 620 enlout(1:18)=enlout(1:18)+enljj(1:18) 621 end do 622 end do 623 ABI_DEALLOCATE(cft) 624 ABI_DEALLOCATE(cfu) 625 ABI_DEALLOCATE(enljj) 626 end if 627 628 end if 629 630 if (paw_opt==3) then 631 632 ! ============== Accumulate contribution to <c|S|c> =============== 633 if (choice==1) then 634 do ispinor=1,nspinor 635 do ia=1,nincat 636 do ilmn=1,nlmn 637 do iplex=1,cplex 638 enlout(1)=enlout(1)+gxfac_sij(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor) 639 end do 640 end do 641 end do 642 end do 643 end if 644 645 ! ============== Accumulate contribution to <c|dS/d_atm_pos|c> =============== 646 if (choice==2.or.choice==23) then 647 ishift=0;if (choice==23) ishift=6 648 do ispinor=1,nspinor 649 do ia=1,nincat 650 enlj(1:3)=zero 651 iashift=3*(ia+ia3-2) 652 do ilmn=1,nlmn 653 do mu=1,3 654 do iplex=1,cplex 655 enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor) 656 end do 657 end do 658 end do 659 enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3) 660 end do 661 end do 662 end if 663 664 ! ============== Accumulate contribution to <c|dS/d_strain|c> =============== 665 if (choice==3.or.choice==23) then 666 enlj(1:6)=zero 667 do ispinor=1,nspinor 668 do ia=1,nincat 669 do ilmn=1,nlmn 670 gxfacj(1:cplex)=gxfac_sij(1:cplex,ilmn,ia,ispinor) 671 do iplex=1,cplex 672 enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor) 673 end do 674 do mu=1,6 675 do iplex=1,cplex 676 enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor) 677 end do 678 end do 679 end do 680 end do 681 end do 682 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 683 end if 684 685 ! ======== Accumulate the contributions of derivatives of <c|S|c> wrt to k ========== 686 if (choice==5) then 687 enlj(1:3)=zero 688 ! If cplex=1, gxfac is real and dgxdt is pure imaginary; thus there is no contribution 689 if(cplex==2)then 690 do ispinor=1,nspinor 691 do ia=1,nincat 692 do ilmn=1,nlmn 693 do mu=1,3 694 do iplex=1,cplex 695 enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor) 696 end do 697 end do 698 end do 699 end do 700 end do 701 end if 702 enlout(1:3)=enlout(1:3)+two*enlj(1:3) 703 end if 704 705 ! ====== Accumulate the contributions of left or right derivatives of <c|S|c> wrt to k ========== 706 ! Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k 707 if (choice==51.or.choice==52) then 708 enlj(1:6)=zero 709 do ispinor=1,nspinor 710 do ia=1,nincat 711 if(cplex==2)then 712 do ilmn=1,nlmn 713 do mu=1,3 714 enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) & 715 & +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) 716 enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) & 717 & -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 718 end do 719 end do 720 else if (cplex_fac==2) then 721 do ilmn=1,nlmn 722 do mu=1,3 723 enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 724 enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 725 end do 726 end do 727 else if (cplex_fac==1) then 728 do ilmn=1,nlmn 729 do mu=1,3 730 enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 731 end do 732 end do 733 end if 734 end do 735 end do 736 if (choice==52) then 737 enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6) 738 end if 739 enlout(1:6)=enlout(1:6)+enlj(1:6) 740 end if 741 742 ! ====== Accumulate contribution to <c|d2S/d_atm_pos d_left_k|c> ========= 743 if (choice==54) then 744 ABI_ALLOCATE(enljj,(18)) 745 do ispinor=1,nspinor 746 do ia=1,nincat 747 enljj(1:18)=zero 748 iashift=18*(ia+ia3-2) 749 if(cplex==2) then 750 do ilmn=1,nlmn 751 mu=1;nu=1 752 do mua=1,3 ! atm. pos 753 do mub=1,3 ! k 754 enljj(nu)=enljj(nu) & 755 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) & 756 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) & 757 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) & 758 & +gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) 759 760 enljj(nu+1)=enljj(nu+1) & 761 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) & 762 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) & 763 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) & 764 & -gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 765 766 mu=mu+1;nu=nu+2 767 end do 768 end do 769 end do 770 ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k 771 else 772 do ilmn=1,nlmn 773 mu=1;nu=1 774 do mua=1,3 ! atm. pos 775 do mub=1,3 ! k 776 enljj(nu+1)=enljj(nu+1) & 777 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) & 778 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 779 mu=mu+1;nu=nu+2 780 end do 781 end do 782 end do 783 end if 784 enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18) 785 end do 786 end do 787 ABI_DEALLOCATE(enljj) 788 end if 789 790 ! ====== Accumulate contribution to <c|d2S/d_dstrain d_right_k|c> ========= 791 if (choice==55) then 792 ABI_ALLOCATE(enljj,(36)) 793 do ispinor=1,nspinor 794 do ia=1,nincat 795 enljj(1:36)=zero;enlj(:)=zero 796 ! If cplex=1, dgxdt is real for strain, pure imaginary for k; 797 ! If cplex_fac=1, dgxdtfac is pure imaginary for k; 798 if(cplex==2.and.cplex_fac==2) then 799 do ilmn=1,nlmn 800 ! First compute 2nd-derivative contribution 801 mu=1 802 do mua=1,6 ! strain (lambda,nu) 803 mua1=alpha(mua) ! (nu) 804 mua2=beta(mua) ! (lambda) 805 do mub=1,3 ! k (mu) 806 muu=3*(gamma(mua1,mub)-1)+mua2 807 mut=3*(gamma(mua2,mub)-1)+mua1 808 d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) & 809 & +d2gxdt(1:cplex,mut,ilmn,ia,ispinor)) 810 enljj(mu)=enljj(mu) & 811 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) & 812 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) & 813 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1)+gxfac_sij(2,ilmn,ia,ispinor)*d2gx(2) 814 enljj(mu+1)=enljj(mu+1) & 815 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) & 816 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) & 817 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(2)-gxfac_sij(2,ilmn,ia,ispinor)*d2gx(1) 818 mu=mu+2 819 end do 820 end do 821 ! Then store 1st-derivative contribution 822 mu=1 823 do nu=1,3 824 enlj(mu )=enlj(mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) & 825 & +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) 826 enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) & 827 & -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 828 mu=mu+2 829 end do 830 end do 831 ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k 832 else 833 do ilmn=1,nlmn 834 mu=1 835 do mua=1,6 ! strain (lambda,nu) 836 mua1=alpha(mua) ! (nu) 837 mua2=beta(mua) ! (lambda) 838 do mub=1,3 ! k (mu) 839 muu=3*(gamma(mua1,mub)-1)+mua2 840 mut=3*(gamma(mua2,mub)-1)+mua1 841 d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor)) 842 enljj(mu+1)=enljj(mu+1) & 843 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) & 844 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1) 845 mu=mu+2 846 end do 847 end do 848 ! Then store 1st-derivative contribution 849 mu=1 850 do nu=1,3 851 enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 852 mu=mu+2 853 end do 854 end do 855 end if 856 enlout(1:36)=enlout(1:36)+enljj(1:36) 857 ddkk(1:6)=ddkk(1:6)+enlj(1:6) 858 end do 859 end do 860 ABI_DEALLOCATE(enljj) 861 end if 862 863 ! ====== Accumulate contribution to <c|d2S/d_k d_k|c> ========= 864 if (choice==8) then 865 ABI_ALLOCATE(cft,(3,nlmn)) 866 ABI_ALLOCATE(cfu,(3,nlmn)) 867 do ispinor=1,nspinor 868 do ia=1,nincat 869 enlj(1:6)=zero 870 do ilmn=1,nlmn 871 do mu=1,6 872 do iplex=1,cplex 873 enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor) 874 end do 875 end do 876 end do 877 ! If cplex=1, dgxdt is pure imaginary, dgxdtfac_sij is pure imaginary; 878 if(cplex==2)then 879 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 880 cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor)) 881 else 882 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 883 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor)) 884 end if 885 do ilmn=1,nlmn 886 do mu=1,6 887 mua=alpha(mu);mub=beta(mu) 888 enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn)) 889 end do 890 end do 891 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 892 end do 893 end do 894 ABI_DEALLOCATE(cft) 895 ABI_DEALLOCATE(cfu) 896 end if 897 898 ! ====== Accumulate contribution to <c|d/d_k[d(right)S/d_k]|c> ========= 899 ! Full derivative wrt to k1, right derivative wrt to k2 900 if (choice==81) then 901 ABI_ALLOCATE(cft,(3,nlmn)) 902 ABI_ALLOCATE(cfu,(6,nlmn)) 903 ABI_ALLOCATE(enljj,(18)) 904 do ispinor=1,nspinor 905 do ia=1,nincat 906 enljj(1:18)=zero 907 if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real 908 cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),gxfac_sij(2,1:nlmn,ia,ispinor)) 909 else 910 cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),zero) 911 end if 912 if(cplex==2)then !If cplex=1, d2gxdt is pure real 913 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor)) 914 else 915 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero) 916 end if 917 do ilmn=1,nlmn 918 do mu=1,3 919 do nu=1,3 920 muu=3*(mu-1)+nu ; mut=gamma(mu,nu) 921 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 922 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 923 end do 924 end do 925 end do 926 if(cplex==2)then !If cplex=1, dgxdt is pure imaginary 927 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 928 else 929 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 930 end if 931 if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary 932 cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor)) 933 else 934 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor)) 935 end if 936 do ilmn=1,nlmn 937 do mu=1,3 938 do nu=1,3 939 muu=3*(mu-1)+nu 940 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 941 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 942 end do 943 end do 944 end do 945 enlout(1:18)=enlout(1:18)+enljj(1:18) 946 end do 947 end do 948 ABI_DEALLOCATE(cft) 949 ABI_DEALLOCATE(cfu) 950 ABI_DEALLOCATE(enljj) 951 end if 952 953 end if 954 955 end subroutine opernld_ylm