TABLE OF CONTENTS
ABINIT/m_opernld_ylm [ Modules ]
NAME
m_opernld_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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_opernld_ylm 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 28 implicit none 29 30 private
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
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=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=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.
SOURCE
120 subroutine opernld_ylm(choice,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,& 121 & enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,ndat_left,nd2gxdt,ndgxdt,& 122 & ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt,strnlk,& 123 & enlout_im) 124 125 !Arguments ------------------------------------ 126 !scalars 127 integer,intent(in) :: choice,cplex,cplex_fac,ia3,natom,nd2gxdt,ndgxdt 128 integer,intent(in) :: ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt 129 integer,intent(in) :: ndat_left 130 real(dp),intent(inout) :: enlk 131 !arrays 132 real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor) 133 real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor) 134 real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor) 135 real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3)) 136 real(dp),intent(in) :: gx(cplex,nlmn,nincat,nspinor*ndat_left),gxfac(cplex_fac,nlmn,nincat,nspinor) 137 real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3)) 138 real(dp),intent(inout) :: ddkk(6),enlout(nnlout*ndat_left),fnlk(3*natom),strnlk(6) 139 real(dp),intent(inout),optional :: enlout_im(nnlout*ndat_left) 140 141 !Local variables------------------------------- 142 !scalars 143 integer :: ia,iashift,idat_left,ilmn,iplex,ishift,ispinor,mu,mua,mua1,mua2,mub,mushift,mut,muu 144 integer :: nu,nushift 145 real(dp) :: dummy 146 !arrays 147 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 148 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 149 integer,parameter :: twist_dir(6)=(/2,3,3,1,1,2/) 150 real(dp) :: d2gx(cplex),enlj(6*cplex),gxfacj(cplex) 151 real(dp),allocatable :: enljj(:) 152 complex(dpc),allocatable :: cft(:,:), cfu(:,:) 153 154 ! ************************************************************************* 155 156 ABI_CHECK(cplex_fac>=cplex,'BUG: invalid cplex_fac<cplex!') 157 158 if (paw_opt==0.or.paw_opt==1.or.paw_opt==2) then 159 160 ! ============== Accumulate the non-local energy =============== 161 if (choice==1) then 162 if (present(enlout_im).and.cplex==2) then ! cplex=cplex_fac=2 163 do idat_left=1,ndat_left 164 do ispinor=1,nspinor 165 do ia=1,nincat 166 do ilmn=1,nlmn 167 enlout (idat_left)=enlout (idat_left)+gxfac(1,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor) 168 enlout (idat_left)=enlout (idat_left)+gxfac(2,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor) 169 enlout_im(idat_left)=enlout_im(idat_left)+gxfac(2,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor) 170 enlout_im(idat_left)=enlout_im(idat_left)-gxfac(1,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor) 171 end do 172 end do 173 end do 174 end do 175 else if (present(enlout_im).and.cplex_fac==2) then ! cplex=1,cplex_fac=2 176 do idat_left=1,ndat_left 177 do ispinor=1,nspinor 178 do ia=1,nincat 179 do ilmn=1,nlmn 180 enlout (idat_left)=enlout (idat_left)+gxfac(1,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor) 181 enlout_im(idat_left)=enlout_im(idat_left)+gxfac(2,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor) 182 end do 183 end do 184 end do 185 end do 186 else ! only the real part is needed or the imaginary part is zero 187 do idat_left=1,ndat_left 188 do ispinor=1,nspinor 189 do ia=1,nincat 190 do ilmn=1,nlmn 191 do iplex=1,cplex 192 enlout(idat_left)=enlout(idat_left)+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor+(idat_left-1)*nspinor) 193 end do 194 end do 195 end do 196 end do 197 end do 198 end if 199 end if 200 201 ! ============ Accumulate the forces contributions ============= 202 if (choice==2.or.choice==23.or.choice==24) then 203 ishift=0;if (choice==23) ishift=6 204 do ispinor=1,nspinor 205 do ia=1,nincat 206 enlj(1:3)=zero 207 iashift=3*(ia+ia3-2)+ishift 208 do ilmn=1,nlmn 209 do mu=1,3 210 dummy = zero ! Dummy needed here to get the correct forces with intel -O3 211 do iplex=1,cplex 212 !enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor) 213 dummy=dummy+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor) 214 end do 215 enlj(mu)=enlj(mu)+dummy 216 end do 217 end do 218 enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3) 219 end do 220 end do 221 end if 222 223 ! ======== Accumulate the stress tensor contributions ========== 224 if (choice==3.or.choice==23) then 225 enlj(1:6)=zero 226 do ispinor=1,nspinor 227 do ia=1,nincat 228 do ilmn=1,nlmn 229 gxfacj(1:cplex)=gxfac(1:cplex,ilmn,ia,ispinor) 230 do iplex=1,cplex 231 enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor) 232 end do 233 do mu=1,6 234 do iplex=1,cplex 235 enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor) 236 end do 237 end do 238 end do 239 end do 240 end do 241 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 242 end if 243 244 ! ====== Accumulate the dynamical matrix contributions ========= 245 if (choice==4.or.choice==24) then 246 ishift=0;if (choice==24) ishift=3*natom 247 do ispinor=1,nspinor 248 do ia=1,nincat 249 enlj(1:6)=zero 250 iashift=6*(ia+ia3-2)+ishift 251 do ilmn=1,nlmn 252 do mu=1,6 253 mua=alpha(mu);mub=beta(mu) 254 do iplex=1,cplex 255 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)& 256 & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,mua,ilmn,ia,ispinor) 257 end do 258 end do 259 end do 260 enlout(iashift+1:iashift+6)=enlout(iashift+1:iashift+6)+two*enlj(1:6) 261 end do 262 end do 263 end if 264 265 ! ======== Accumulate the contributions of derivatives of E wrt to k ========== 266 if (choice==5) then 267 enlj(1:3)=zero 268 do ispinor=1,nspinor 269 do ia=1,nincat 270 if(cplex==2)then 271 do ilmn=1,nlmn 272 do mu=1,3 273 do iplex=1,cplex 274 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor) 275 end do 276 end do 277 end do 278 ! If cplex=1, dgxdt is pure imaginary; thus there is no contribution 279 else if (cplex_fac==2) then 280 do ilmn=1,nlmn 281 do mu=1,3 282 enlj(mu)=enlj(mu)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 283 end do 284 end do 285 end if 286 end do 287 end do 288 enlout(1:3)=enlout(1:3)+two*enlj(1:3) 289 end if 290 291 ! ======== Accumulate the contributions of partial derivatives of E wrt to k ========== 292 ! Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k 293 if (choice==51.or.choice==52) then 294 enlj(1:6)=zero 295 do ispinor=1,nspinor 296 do ia=1,nincat 297 if(cplex==2)then 298 do ilmn=1,nlmn 299 do mu=1,3 300 enlj(2*mu-1)=enlj(2*mu-1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) & 301 & +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) 302 enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) & 303 & -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 304 end do 305 end do 306 else if (cplex_fac==2) then 307 do ilmn=1,nlmn 308 do mu=1,3 309 enlj(2*mu-1)=enlj(2*mu-1)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 310 enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 311 end do 312 end do 313 else if (cplex_fac==1) then 314 do ilmn=1,nlmn 315 do mu=1,3 316 enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 317 end do 318 end do 319 end if 320 end do 321 end do 322 if (choice==52) then 323 enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6) 324 end if 325 enlout(1:6)=enlout(1:6)+enlj(1:6) 326 end if 327 328 ! ======== Accumulate the contributions of twist derivatives of E wrt to k ========== 329 ! accumulate <u|dp_i/dk_(idir+1)>D_ij<dp_j/dk(idir+2)|u> 330 if (choice==53) then 331 enlj(:)=zero 332 ABI_MALLOC(cft,(3,nlmn)) 333 ABI_MALLOC(cfu,(3,nlmn)) 334 ! If cplex=1, dgxdt is pure imaginary; 335 ! If cplex_fac=1, dgxdtfac is pure imaginary; 336 do ispinor=1,nspinor 337 do ia=1,nincat 338 if(cplex==2)then 339 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 340 else 341 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 342 end if 343 if(cplex_fac==2)then 344 cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor)) 345 else 346 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor)) 347 end if 348 do ilmn=1,nlmn 349 do mu=1,3 350 mut = twist_dir(2*mu-1) 351 muu = twist_dir(2*mu) 352 if (cplex == 2) then 353 enlj(2*mu-1) = enlj(2*mu-1) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 354 enlj(2*mu) = enlj(2*mu) + aimag(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 355 else 356 enlj(mu) = enlj(mu) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 357 end if 358 end do ! end loop over mu=1,3 359 end do ! end loop over ilmn states 360 end do ! end loop over ia atoms 361 end do ! end loop over ispinor 362 do mu = 1, 3 363 if (cplex == 2) then 364 enlout(2*mu-1)=enlout(2*mu-1)+enlj(2*mu-1) 365 enlout(2*mu) =enlout(2*mu) +enlj(2*mu) 366 else 367 enlout(mu)=enlout(mu)+enlj(mu) 368 end if 369 end do ! end loop over mu = 1, 3 370 ABI_FREE(cft) 371 ABI_FREE(cfu) 372 end if 373 374 ! ====== Accumulate the effective charges contributions ========= 375 if (choice==54) then 376 ABI_MALLOC(enljj,(18)) 377 do ispinor=1,nspinor 378 do ia=1,nincat 379 enljj(1:18)=zero 380 iashift=18*(ia+ia3-2) 381 ! If cplex=1, dgxdt is real for atm. pos, pure imaginary for k; 382 ! If cplex_fac=1, dgxdtfac is pure imaginary for k; 383 if(cplex==2.and.cplex_fac==2) then 384 do ilmn=1,nlmn 385 mu=1;nu=1 386 do mua=1,3 ! atm. pos 387 do mub=1,3 ! k 388 enljj(nu)=enljj(nu) & 389 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 390 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) & 391 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) & 392 & +gxfac(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) 393 enljj(nu+1)=enljj(nu+1) & 394 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) & 395 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 396 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) & 397 & -gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 398 mu=mu+1;nu=nu+2 399 end do 400 end do 401 end do 402 else if(cplex==1.and.cplex_fac==2)then 403 do ilmn=1,nlmn 404 mu=1;nu=1 405 do mua=1,3 ! atm. pos 406 do mub=1,3 ! k 407 enljj(nu)=enljj(nu) & 408 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 409 & +gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 410 enljj(nu+1)=enljj(nu+1) & 411 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) & 412 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 413 mu=mu+1;nu=nu+2 414 end do 415 end do 416 end do 417 else if(cplex==1.and.cplex_fac==1)then 418 do ilmn=1,nlmn 419 mu=1;nu=1 420 do mua=1,3 ! atm. pos 421 do mub=1,3 ! k 422 enljj(nu+1)=enljj(nu+1) & 423 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) & 424 & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 425 mu=mu+1;nu=nu+2 426 end do 427 end do 428 end do 429 end if 430 enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18) 431 end do 432 end do 433 ABI_FREE(enljj) 434 end if 435 436 ! ====== Accumulate the piezoelectric tensor contributions ========= 437 if (choice==55) then 438 ABI_MALLOC(enljj,(36)) 439 do ispinor=1,nspinor 440 do ia=1,nincat 441 enljj(1:36)=zero;enlj(:)=zero 442 ! If cplex=1, dgxdt is real for strain, pure imaginary for k; 443 ! If cplex_fac=1, dgxdtfac is pure imaginary for k; 444 if(cplex==2.and.cplex_fac==2) then 445 do ilmn=1,nlmn 446 ! First compute 2nd-derivative contribution 447 mu=1 448 do mua=1,6 ! strain (lambda,nu) 449 mua1=alpha(mua) ! (nu) 450 mua2=beta(mua) ! (lambda) 451 do mub=1,3 ! k (mu) 452 muu=3*(gamma(mua1,mub)-1)+mua2 453 mut=3*(gamma(mua2,mub)-1)+mua1 454 d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) & 455 & +d2gxdt(1:cplex,mut,ilmn,ia,ispinor)) 456 enljj(mu)=enljj(mu) & 457 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 458 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) & 459 & +gxfac(1,ilmn,ia,ispinor)*d2gx(1)+gxfac(2,ilmn,ia,ispinor)*d2gx(2) 460 enljj(mu+1)=enljj(mu+1) & 461 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) & 462 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 463 & +gxfac(1,ilmn,ia,ispinor)*d2gx(2)-gxfac(2,ilmn,ia,ispinor)*d2gx(1) 464 mu=mu+2 465 end do 466 end do 467 ! Then store 1st-derivative contribution 468 mu=1 469 do nu=1,3 470 enlj(mu )=enlj(mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) & 471 & +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) 472 enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) & 473 & -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 474 mu=mu+2 475 end do 476 end do 477 else if(cplex==1.and.cplex_fac==2)then 478 do ilmn=1,nlmn 479 ! First compute 2nd-derivative contribution 480 mu=1 481 do mua=1,6 ! strain (lambda,nu) 482 mua1=alpha(mua) ! (nu) 483 mua2=beta(mua) ! (lambda) 484 do mub=1,3 ! k (mu) 485 muu=3*(gamma(mua1,mub)-1)+mua2 486 mut=3*(gamma(mua2,mub)-1)+mua1 487 d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor)) 488 enljj(mu)=enljj(mu) & 489 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 490 & +gxfac(2,ilmn,ia,ispinor)*d2gx(1) 491 enljj(mu+1)=enljj(mu+1) & 492 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) & 493 & +gxfac(1,ilmn,ia,ispinor)*d2gx(1) 494 mu=mu+2 495 end do 496 end do 497 ! Then store 1st-derivative contribution 498 mu=1 499 do nu=1,3 500 enlj(mu )=enlj(mu )+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 501 enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 502 mu=mu+2 503 end do 504 end do 505 else if(cplex==1.and.cplex_fac==1)then 506 do ilmn=1,nlmn 507 mu=1 508 do mua=1,6 ! strain (lambda,nu) 509 mua1=alpha(mua) ! (nu) 510 mua2=beta(mua) ! (lambda) 511 do mub=1,3 ! k (mu) 512 muu=3*(gamma(mua1,mub)-1)+mua2 513 mut=3*(gamma(mua2,mub)-1)+mua1 514 d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor)) 515 enljj(mu+1)=enljj(mu+1) & 516 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) & 517 & +gxfac(1,ilmn,ia,ispinor)*d2gx(1) 518 mu=mu+2 519 end do 520 end do 521 ! Then store 1st-derivative contribution 522 mu=1 523 do nu=1,3 524 enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 525 mu=mu+2 526 end do 527 end do 528 end if 529 enlout(1:36)=enlout(1:36)+enljj(1:36) 530 ddkk(1:6)=ddkk(1:6)+enlj(1:6) 531 end do 532 end do 533 ABI_FREE(enljj) 534 end if 535 536 ! ======= Accumulate the elastic tensor contributions ========== 537 if (choice==6) then 538 do ispinor=1,nspinor 539 do ia=1,nincat 540 iashift=3*(ia+ia3-2) 541 do ilmn=1,nlmn 542 do iplex=1,cplex 543 enlk=enlk+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor) 544 end do 545 enlj(1:3)=zero 546 do mu=1,3 547 do iplex=1,cplex 548 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,6+mu,ilmn,ia,ispinor) 549 end do 550 end do 551 fnlk(iashift+1:iashift+3)=fnlk(iashift+1:iashift+3)+two*enlj(1:3) 552 enlj(1:6)=zero 553 do mu=1,6 554 do iplex=1,cplex 555 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor) 556 end do 557 end do 558 strnlk(1:6)=strnlk(1:6)+two*enlj(1:6) 559 do mub=1,6 560 mushift=6*(mub-1);nushift=(3*natom+6)*(mub-1) 561 do mua=1,6 562 mu=mushift+mua;nu=nushift+mua 563 do iplex=1,cplex 564 enlout(nu)=enlout(nu)+two* & 565 & (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)& 566 & +dgxdtfac(iplex,mua,ilmn,ia,ispinor)*dgxdt(iplex,mub,ilmn,ia,ispinor)) 567 end do 568 end do 569 mushift=36+3*(mub-1);nushift=6+iashift+(3*natom+6)*(mub-1) 570 do mua=1,3 571 mu=mushift+mua;nu=nushift+mua 572 do iplex=1,cplex 573 enlout(nu)=enlout(nu)+two* & 574 & (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)& 575 & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,6+mua,ilmn,ia,ispinor)) 576 end do 577 end do 578 end do 579 end do 580 end do 581 end do 582 end if 583 584 ! ======== Accumulate the contributions of 2nd-derivatives of E wrt to k ========== 585 if (choice==8) then 586 ABI_MALLOC(cft,(3,nlmn)) 587 ABI_MALLOC(cfu,(3,nlmn)) 588 do ispinor=1,nspinor 589 do ia=1,nincat 590 enlj(1:6)=zero 591 do ilmn=1,nlmn 592 do mu=1,6 593 do iplex=1,cplex 594 enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor) 595 end do 596 end do 597 end do 598 ! If cplex=1, dgxdt is pure imaginary; 599 ! If cplex_fac=1, dgxdtfac is pure imaginary; 600 if(cplex==2)then 601 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 602 else 603 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 604 end if 605 if(cplex_fac==2)then 606 cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor)) 607 else 608 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor)) 609 end if 610 do ilmn=1,nlmn 611 do mu=1,6 612 mua=alpha(mu);mub=beta(mu) 613 enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn)) 614 end do 615 end do 616 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 617 end do 618 end do 619 ABI_FREE(cft) 620 ABI_FREE(cfu) 621 end if 622 623 ! ======== Accumulate the contributions of partial 2nd-derivatives of E wrt to k ========== 624 ! Full derivative wrt to k1, right derivative wrt to k2 625 if (choice==81) then 626 ABI_MALLOC(cft,(3,nlmn)) 627 ABI_MALLOC(cfu,(6,nlmn)) 628 ABI_MALLOC(enljj,(18)) 629 do ispinor=1,nspinor 630 do ia=1,nincat 631 enljj(1:18)=zero 632 if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real 633 cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),gxfac(2,1:nlmn,ia,ispinor)) 634 else 635 cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),zero) 636 end if 637 if(cplex==2)then !If cplex=1, d2gxdt is pure real 638 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor)) 639 else 640 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero) 641 end if 642 do ilmn=1,nlmn 643 do mu=1,3 644 do nu=1,3 645 muu=3*(mu-1)+nu ; mut=gamma(mu,nu) 646 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 647 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 648 end do 649 end do 650 end do 651 if(cplex==2)then !If cplex=1, dgxdt is pure imaginary 652 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 653 else 654 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 655 end if 656 if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary 657 cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor)) 658 else 659 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor)) 660 end if 661 do ilmn=1,nlmn 662 do mu=1,3 663 do nu=1,3 664 muu=3*(mu-1)+nu 665 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 666 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 667 end do 668 end do 669 end do 670 enlout(1:18)=enlout(1:18)+enljj(1:18) 671 end do 672 end do 673 ABI_FREE(cft) 674 ABI_FREE(cfu) 675 ABI_FREE(enljj) 676 end if 677 678 end if 679 680 if (paw_opt==3) then 681 682 ! ============== Accumulate contribution to <c|S|c> =============== 683 if (choice==1) then 684 if (present(enlout_im).and.cplex==2) then ! cplex=2 685 do idat_left=1,ndat_left 686 do ispinor=1,nspinor 687 do ia=1,nincat 688 do ilmn=1,nlmn 689 enlout (idat_left)=enlout (idat_left)+gxfac_sij(1,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor) 690 enlout (idat_left)=enlout (idat_left)+gxfac_sij(2,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor) 691 enlout_im(idat_left)=enlout_im(idat_left)+gxfac_sij(2,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor) 692 enlout_im(idat_left)=enlout_im(idat_left)-gxfac_sij(1,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor) 693 end do 694 end do 695 end do 696 end do 697 else ! only the real part is needed or the imaginary part is zero 698 do idat_left=1,ndat_left 699 do ispinor=1,nspinor 700 do ia=1,nincat 701 do ilmn=1,nlmn 702 do iplex=1,cplex 703 enlout(idat_left)=enlout(idat_left)+& 704 & gxfac_sij(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor+(idat_left-1)*nspinor) 705 end do 706 end do 707 end do 708 end do 709 end do 710 end if 711 end if 712 713 ! ============== Accumulate contribution to <c|dS/d_atm_pos|c> =============== 714 if (choice==2.or.choice==23) then 715 ishift=0;if (choice==23) ishift=6 716 do ispinor=1,nspinor 717 do ia=1,nincat 718 enlj(1:3)=zero 719 iashift=3*(ia+ia3-2) 720 do ilmn=1,nlmn 721 do mu=1,3 722 do iplex=1,cplex 723 enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor) 724 end do 725 end do 726 end do 727 enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3) 728 end do 729 end do 730 end if 731 732 ! ============== Accumulate contribution to <c|dS/d_strain|c> =============== 733 if (choice==3.or.choice==23) then 734 enlj(1:6)=zero 735 do ispinor=1,nspinor 736 do ia=1,nincat 737 do ilmn=1,nlmn 738 gxfacj(1:cplex)=gxfac_sij(1:cplex,ilmn,ia,ispinor) 739 do iplex=1,cplex 740 enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor) 741 end do 742 do mu=1,6 743 do iplex=1,cplex 744 enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor) 745 end do 746 end do 747 end do 748 end do 749 end do 750 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 751 end if 752 753 ! ======== Accumulate the contributions of derivatives of <c|S|c> wrt to k ========== 754 if (choice==5) then 755 enlj(1:3)=zero 756 ! If cplex=1, gxfac is real and dgxdt is pure imaginary; thus there is no contribution 757 if(cplex==2)then 758 do ispinor=1,nspinor 759 do ia=1,nincat 760 do ilmn=1,nlmn 761 do mu=1,3 762 do iplex=1,cplex 763 enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor) 764 end do 765 end do 766 end do 767 end do 768 end do 769 end if 770 enlout(1:3)=enlout(1:3)+two*enlj(1:3) 771 end if 772 773 ! ====== Accumulate the contributions of left or right derivatives of <c|S|c> wrt to k ========== 774 ! Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k 775 if (choice==51.or.choice==52) then 776 enlj(1:6)=zero 777 do ispinor=1,nspinor 778 do ia=1,nincat 779 if(cplex==2)then 780 do ilmn=1,nlmn 781 do mu=1,3 782 enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) & 783 & +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) 784 enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) & 785 & -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 786 end do 787 end do 788 else if (cplex_fac==2) then 789 do ilmn=1,nlmn 790 do mu=1,3 791 enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 792 enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 793 end do 794 end do 795 else if (cplex_fac==1) then 796 do ilmn=1,nlmn 797 do mu=1,3 798 enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) 799 end do 800 end do 801 end if 802 end do 803 end do 804 if (choice==52) then 805 enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6) 806 end if 807 enlout(1:6)=enlout(1:6)+enlj(1:6) 808 end if 809 810 ! ====== Accumulate the contributions of twist derivatives of <c|S|c> wrt to k ========== 811 ! Choice 53: <u|dp_i/dk_(idir+1)>Sij<dp_j/dk_(idir+2)|u> 812 if (choice==53) then 813 enlj(:)=zero 814 ABI_MALLOC(cft,(3,nlmn)) 815 ABI_MALLOC(cfu,(3,nlmn)) 816 ! If cplex=1, dgxdt is pure imaginary; 817 ! If cplex_fac=1, dgxdtfac is pure imaginary; 818 do ispinor=1,nspinor 819 do ia=1,nincat 820 if(cplex==2)then 821 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 822 else 823 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 824 end if 825 if(cplex_fac==2)then 826 cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor)) 827 else 828 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor)) 829 end if 830 do ilmn=1,nlmn 831 do mu=1,3 832 mut = twist_dir(2*mu-1) 833 muu = twist_dir(2*mu) 834 if (cplex == 2) then 835 enlj(2*mu-1) = enlj(2*mu-1) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 836 enlj(2*mu) = enlj(2*mu) + aimag(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 837 else 838 enlj(mu) = enlj(mu) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn)) 839 end if 840 end do ! end loop over mu=1,3 841 end do ! end loop over ilmn states 842 end do ! end loop over ia atoms 843 end do ! end loop over ispinor 844 do mu = 1, 3 845 if (cplex == 2) then 846 enlout(2*mu-1)=enlout(2*mu-1)+enlj(2*mu-1) 847 enlout(2*mu) =enlout(2*mu) +enlj(2*mu) 848 else 849 enlout(mu)=enlout(mu)+enlj(mu) 850 end if 851 end do ! end loop over mu = 1, 3 852 ABI_FREE(cft) 853 ABI_FREE(cfu) 854 end if 855 856 ! ====== Accumulate contribution to <c|d2S/d_atm_pos d_left_k|c> ========= 857 if (choice==54) then 858 ABI_MALLOC(enljj,(18)) 859 do ispinor=1,nspinor 860 do ia=1,nincat 861 enljj(1:18)=zero 862 iashift=18*(ia+ia3-2) 863 if(cplex==2) then 864 do ilmn=1,nlmn 865 mu=1;nu=1 866 do mua=1,3 ! atm. pos 867 do mub=1,3 ! k 868 enljj(nu)=enljj(nu) & 869 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) & 870 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) & 871 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) & 872 & +gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) 873 874 enljj(nu+1)=enljj(nu+1) & 875 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) & 876 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) & 877 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) & 878 & -gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 879 880 mu=mu+1;nu=nu+2 881 end do 882 end do 883 end do 884 ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k 885 else 886 do ilmn=1,nlmn 887 mu=1;nu=1 888 do mua=1,3 ! atm. pos 889 do mub=1,3 ! k 890 enljj(nu+1)=enljj(nu+1) & 891 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) & 892 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) 893 mu=mu+1;nu=nu+2 894 end do 895 end do 896 end do 897 end if 898 enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18) 899 end do 900 end do 901 ABI_FREE(enljj) 902 end if 903 904 ! ====== Accumulate contribution to <c|d2S/d_dstrain d_right_k|c> ========= 905 if (choice==55) then 906 ABI_MALLOC(enljj,(36)) 907 do ispinor=1,nspinor 908 do ia=1,nincat 909 enljj(1:36)=zero;enlj(:)=zero 910 ! If cplex=1, dgxdt is real for strain, pure imaginary for k; 911 ! If cplex_fac=1, dgxdtfac is pure imaginary for k; 912 if(cplex==2.and.cplex_fac==2) then 913 do ilmn=1,nlmn 914 ! First compute 2nd-derivative contribution 915 mu=1 916 do mua=1,6 ! strain (lambda,nu) 917 mua1=alpha(mua) ! (nu) 918 mua2=beta(mua) ! (lambda) 919 do mub=1,3 ! k (mu) 920 muu=3*(gamma(mua1,mub)-1)+mua2 921 mut=3*(gamma(mua2,mub)-1)+mua1 922 d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) & 923 & +d2gxdt(1:cplex,mut,ilmn,ia,ispinor)) 924 enljj(mu)=enljj(mu) & 925 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) & 926 & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) & 927 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1)+gxfac_sij(2,ilmn,ia,ispinor)*d2gx(2) 928 enljj(mu+1)=enljj(mu+1) & 929 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) & 930 & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) & 931 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(2)-gxfac_sij(2,ilmn,ia,ispinor)*d2gx(1) 932 mu=mu+2 933 end do 934 end do 935 ! Then store 1st-derivative contribution 936 mu=1 937 do nu=1,3 938 enlj(mu )=enlj(mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) & 939 & +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) 940 enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) & 941 & -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 942 mu=mu+2 943 end do 944 end do 945 ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k 946 else 947 do ilmn=1,nlmn 948 mu=1 949 do mua=1,6 ! strain (lambda,nu) 950 mua1=alpha(mua) ! (nu) 951 mua2=beta(mua) ! (lambda) 952 do mub=1,3 ! k (mu) 953 muu=3*(gamma(mua1,mub)-1)+mua2 954 mut=3*(gamma(mua2,mub)-1)+mua1 955 d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor)) 956 enljj(mu+1)=enljj(mu+1) & 957 & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) & 958 & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1) 959 mu=mu+2 960 end do 961 end do 962 ! Then store 1st-derivative contribution 963 mu=1 964 do nu=1,3 965 enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) 966 mu=mu+2 967 end do 968 end do 969 end if 970 enlout(1:36)=enlout(1:36)+enljj(1:36) 971 ddkk(1:6)=ddkk(1:6)+enlj(1:6) 972 end do 973 end do 974 ABI_FREE(enljj) 975 end if 976 977 ! ====== Accumulate contribution to <c|d2S/d_k d_k|c> ========= 978 if (choice==8) then 979 ABI_MALLOC(cft,(3,nlmn)) 980 ABI_MALLOC(cfu,(3,nlmn)) 981 do ispinor=1,nspinor 982 do ia=1,nincat 983 enlj(1:6)=zero 984 do ilmn=1,nlmn 985 do mu=1,6 986 do iplex=1,cplex 987 enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor) 988 end do 989 end do 990 end do 991 ! If cplex=1, dgxdt is pure imaginary, dgxdtfac_sij is pure imaginary; 992 if(cplex==2)then 993 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 994 cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor)) 995 else 996 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 997 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor)) 998 end if 999 do ilmn=1,nlmn 1000 do mu=1,6 1001 mua=alpha(mu);mub=beta(mu) 1002 enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn)) 1003 end do 1004 end do 1005 enlout(1:6)=enlout(1:6)+two*enlj(1:6) 1006 end do 1007 end do 1008 ABI_FREE(cft) 1009 ABI_FREE(cfu) 1010 end if 1011 1012 ! ====== Accumulate contribution to <c|d/d_k[d(right)S/d_k]|c> ========= 1013 ! Full derivative wrt to k1, right derivative wrt to k2 1014 if (choice==81) then 1015 ABI_MALLOC(cft,(3,nlmn)) 1016 ABI_MALLOC(cfu,(6,nlmn)) 1017 ABI_MALLOC(enljj,(18)) 1018 do ispinor=1,nspinor 1019 do ia=1,nincat 1020 enljj(1:18)=zero 1021 if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real 1022 cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),gxfac_sij(2,1:nlmn,ia,ispinor)) 1023 else 1024 cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),zero) 1025 end if 1026 if(cplex==2)then !If cplex=1, d2gxdt is pure real 1027 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor)) 1028 else 1029 cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero) 1030 end if 1031 do ilmn=1,nlmn 1032 do mu=1,3 1033 do nu=1,3 1034 muu=3*(mu-1)+nu ; mut=gamma(mu,nu) 1035 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 1036 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn)) 1037 end do 1038 end do 1039 end do 1040 if(cplex==2)then !If cplex=1, dgxdt is pure imaginary 1041 cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor)) 1042 else 1043 cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor)) 1044 end if 1045 if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary 1046 cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor)) 1047 else 1048 cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor)) 1049 end if 1050 do ilmn=1,nlmn 1051 do mu=1,3 1052 do nu=1,3 1053 muu=3*(mu-1)+nu 1054 enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 1055 enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn)) 1056 end do 1057 end do 1058 end do 1059 enlout(1:18)=enlout(1:18)+enljj(1:18) 1060 end do 1061 end do 1062 ABI_FREE(cft) 1063 ABI_FREE(cfu) 1064 ABI_FREE(enljj) 1065 end if 1066 1067 end if 1068 1069 end subroutine opernld_ylm